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ABSTRACT 


Direct-Sequence Spread-Spectrum (DS-SS) is among the preferred modulation 
techniques for military applications. DS-SS offers three greatly desired characteristics. It 
allows for the development of Low Probability of Detection (LPD) and Low Probability 
of Intercept (LPI) systems and has a very good performance in fading channels. This the- 
sis investigates the performance of the “Cross-Product RV (CRV) decomposition” as the 
basis of blind-equalization algorithms. The CRV is a rank-revealing decomposition alter- 
native to the Eigenvalue Decomposition (EVD) that can provide a recursively updated 
estimate of the signal and noise subspace at a reduced computational cost. The CRV up- 
dating algorithm is implemented in MATLAB and evaluated in a previously proposed 
communication scheme intended for use in an underwater acoustic network called Sea- 
web. The underwater channel is modeled with the Monterey-Miami Parabolic Equation 
Model (MMPE) for various multipath perturbations. The receiver performance is exam- 
ined using a Monte Carlo simulation. Bit-error rates versus signal-to-noise ratio are pre- 


sented for various, noise assumptions, and receiver synchronization assumptions. 
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EXECUTIVE SUMMARY 


Underwater wireless networks have many military and commercial applications. 
With such networks, the control and coordination of various assets, such as mini- 
submarines, divers, submerged instruments and Unmanned Underwater Vehicles (UUVs) 


is possible without the burden of cables or human intervention. 


The underwater environment is possibly the most challenging environment for 
wireless communications. The propagation of sound is greatly affected by transmission 
loss (TL), ambient noise, reverberation, and the overall variability of the channel. These 
limiting factors impose various trade-offs in the design of a reliable communication sys- 
tem. A higher transmission rate implies greater bandwidth, higher TL and intersymbol 


interference (ISI) but an almost negligible channel variation over a single-bit interval. 


Direct-sequence spread-spectrum (DS-SS) is one of the preferred communication 
techniques in military applications. Spread-spectrum modulation has a very good per- 
formance in multipath environments and allows for the development of Low Probability 
of Detection (LPD) and Low Probability of Intercept (LPI) systems. Because of its de- 


sired characteristics, DS-SS is well suited for shallow-water communications. 


This thesis investigates the performance of DS-SS with blind-equalization in 
simulated ocean channels. The equalization algorithms examined are based on the 
“Cross-Product RV (CRV),” a relatively new subspace decomposition algorithm pro- 
posed as an alternative to the Eigenvalue decomposition (EVD). The CRV updating algo- 
rithm provides a recursively updated estimate of the signal and noise subspace at less 
computational cost than the EVD. The transmitter and receiver structures are imple- 
mented in MATLAB and the channel impulse response is simulated with the Monterey- 
Miami Parabolic Equation model (MMPE). The impulse responses extracted account for 
various channel perturbations (interface roughness, internal waves, etc.) and also for the 


Doppler effect caused by the source motion. 


The performance of the proposed communication scheme is evaluated using 


Monte Carlo simulation for various signal-to-noise ratios (SNR) and design parameters 


XVil 


(packet size, data rate and channel length). For all the individual perturbations simulated 
with the MMPE model a Bit Error Rate (BER) below 10° is achievable for a relatively 
low SNR. In the case where all perturbations are combined or synchronization between 


the receiver and the transmitter is not achieved, the receiver’s performance is limited. 


XViil 


I. INTRODUCTION 


A. BACKGROUND 

The primary purpose of this thesis is to improve upon two previously proposed 
Direct-Sequence Spread-Spectrum (DS-SS) communication schemes [1, 2] intended for 
use in an underwater acoustic communication network known as Seaweb. Seaweb is be- 
ing developed by the Space and Naval Warfare Systems Center, San Diego, [3] and is 
motivated by a requirement for wide-area undersea surveillance in littoral waters. This 
undersea wireless network provides the command, control and communications infra- 
structure [4] for coordinating appropriate assets, such as mini-submarines, divers and 
Unmanned Underwater Vehicles (UUVs) to accomplish a given mission in an arbitrary 
ocean environment. The first underwater communication scheme [1] involved a Direct- 
Sequence Differential Binary Phase-Shift Keying with quadrature spreading (DS-IQ- 
DBSK). Error-correction coding and a RAKE receiver for multipath reception were used 
to improve system performance. The second and latest communication scheme as part of 
a previous thesis research [2] tried to develop an improved receiver structure based on 
blind-equalization algorithms. In both schemes [1, 2] the impulse response of the under- 
water acoustic channel was generated using the Bel/hop acoustic propagation model [5]. 
The Bellhop model is a “variant of Gaussian beam tracing and is especially suited for 
high-frequency deep-water problems where normal mode and parabolic equation (PE) 
models are not practical [6]”. In this thesis, the channel impulse response is estimated us- 
ing the Monterey-Miami Parabolic Equation (MMPE) model described in [7]. This is a 
range dependent wave-theory model based upon the parabolic approximation of the wave 
equation [6]. It is a widely used numerical tool for solving the acoustic wave equation 
especially for low frequencies and/or range-dependent environments where it is computa- 
tionally more efficient. 
B. SEAWEB UNDERWATER NETWORK REQUIREMENTS 

Like all other military applications, Seaweb must ensure communication security 
in two ways. It must limit the ability of unauthorized users to detect the presence of the 
communication signal (i.e., low probability of detection) and at the same time decrease 


the ability of a hostile observer to “listen in” when communications are taking place (i.e., 


1 


low probability of intercept). The aforementioned requirements are met by using DS-SS 
as the modulation scheme for packet transmission. DS spread-spectrum systems are able 
to hide the signal in the background noise by operating at low bit-energy per noise- 


spectral-density (z, / N,). Also due to the pseudo-random properties of the code (chip- 


ping sequence) used to spread the message signal, DS systems have a very low probabil- 
ity of intercept. The message can be demodulated only by the intended receivers who 


know the exact replica of the code used at the transmitter [8]. 


Other Seaweb requirements include utility packets of fixed length (72 bits) and a 
limited operating bandwidth of 5 kHz (9-14 kHz operating frequencies). The desired in- 
formation transmission rate is 40 to 100 bits per second (bps) at a bit-error rate (BER) of 
10°. The communication range is between three to five kilometers and the operating 
depth is on the order of 50 to 200 m. 

C. GOALS AND METHODOLOGY 

The first step in improving the previously proposed signaling schemes is to use a 
more accurate underwater acoustic model to account for the temporal and spatial variabil- 
ity of the channel. As already mentioned, this model is the MMPE model and will be de- 
scribed in more detail in Chapter III of this thesis. For the sake of comparison various 


impulse responses are extracted corresponding to different undersea environments. 


The combination of the DS-SS transmitting scheme and the receiver structure 
based on blind equalization algorithms was proven very robust as part of a previous thesis 
[2]. However, the algorithm worked completely offline and the whole received packet 
was needed in order to calculate the filter coefficients and to demodulate the transmitted 
signal efficiently. As part of this thesis research a relatively new subspace tracking algo- 
rithm, the cross product RV, also known as CRV [9], is used to allow for a real-time es- 
timation of the desired subspaces. The transceiver with the CRV algorithm is imple- 
mented in MATLAB for various combinations of multipath, white and colored noise. 

D. BENEFITS OF STUDY 

Underwater wireless communications are much desired in military applications. 

The possibility to maintain signal transmission without the burden of cables enables the 


operation of underwater robots and vehicles (UUVs) and the gathering of data from sub- 


2 


merged instruments without human intervention. An underwater mobile network such as 
Seaweb has many military and commercial applications (submarine communications, 
ocean exploration, etc.). 
E. THESIS ORGANIZATION 

The remainder of this thesis is organized into five chapters. Chapter II discusses 
the general characteristics of underwater acoustic communication channels. Chapter III 
provides a description of the Monterey-Miami Parabolic Equation (MMPE) model and an 
outline of the theory behind the simulated perturbations. Chapter IV describes the CRV 
algorithm, proposed as a substitute of the Eigenvalue Decomposition (EVD) used in [2]. 
The CRV updating algorithm can provide a recursively updated estimate of the signal and 
noise subspace. Chapter V reviews the previously proposed signaling scheme [2] and 
presents the new scheme based on the CRV algorithm. The performance of the new 
transmitter/receiver structure is evaluated under different channel conditions and design 
parameters. Finally, Chapter VI reviews and summarizes the important results and rec- 


ommends follow-on work. 
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Il. CHARACTERISTICS OF UNDERWATER ACOUSTIC 
COMMUNICATION CHANNELS 


The sea forms a remarkably complex waveguide medium for the propagation of 
sound. The purpose of this chapter is to describe the dominant factors affecting sound 
propagation in shallow water and the difficulties those factors impose on designing a reli- 
able communication system. The observation of sound waves is greatly affected by 
transmission loss (TL), ambient noise, and the temporal and spatial variability of the 
channel. Transmission loss and noise determine the signal-to-noise ratio (SNR), the 
maximum transmission range, and the available communication bandwidth. The temporal 
and spatial variability of the channel induced intersymbol interference (ISI), which influ- 
ences system design (i.e., the choice of the modulation/signal detection method and the 
transmitter/receiver structure). 

A. TRANSMISSION LOSS (TL) 

The weakening of an acoustic signal as it travels through the sea is quantitatively 

described by transmission loss, TL. According to [10] and [11], TL is defined as the ratio 


in decibels between the acoustic intensity /(r) in W/m’ at a field point (of distance r) 

and the intensity /,(7%) measured at a distance 7, = 1m from the sound source (transmit- 
ter), 

TL [GB ref. o1=-1010g| “| 01) 2 (2.1) 

For both plane and spherical waves, the acoustic intensity is proportional to the 


square of the pressure amplitude; therefore, TL can also be expressed as 





TL [dB ref. r,]= -20108|# 2 = 20108| (2.2) 
; r 


Transmission loss can be decomposed into two terms, TL caused by en- 


Spreading ? 


ergy spreading and TL due to sound attenuation. Energy spreading depends on 


Attenuation ? 
the propagation distance and the underwater environment characteristics. For the case of 


an omni-directional sound source in an unbounded medium, sound energy undergoes 


spherical spreading and TL due to spreading is proportional to the square of the transmis- 


sion range, 


pee Isospeed [dB ref. 1 m] = 20 log r. (2.3) 


In isospeed underwater channels where the medium is bounded by two plane and 
parallel surfaces, the TL in the far-field is proportional to the first power of the propaga- 


tion distance (cylindrical spreading), 


TL , (dB ref. 1 m]=10logr. (2.4) 


Far-field, Bounded, Isospee 
In shallow water channels acoustic energy undergoes both types of spreading 

mentioned above-spherical spreading close to the source and cylindrical spreading at 

great distances [12]. Thus a more practical expression for the TL in shallow water is 

given by [10] 

20logr fey 


; (2.5) 


TL 
10logr+10logr r>r 


Shallow Water Spreading 


[dB ref. 1 m] -| 


The transition range 7, (i.e., the range at which the change from spherical to cylindrical 


spreading occurs) for the case of an isospeed shallow water channel with fast bottom, is 


estimated by [10] 
r=—, (2.6) 


where #7 is the bottom depth in meters and @. is the critical angle in degrees. 


Shallow Water Channel 











Spherical Spreading Cylindrical Spreading 


Figure 1. Spreading of Acoustic Energy in a Shallow Water Channel. 
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Attenuation of the field is the effect of sound absorption. Absorption involves the 
process of converting acoustic energy into heat, thus resulting in a true loss of energy. 
The dominant causes of absorption in seawater are the ionic relaxation of magnesium sul- 
fate (MgSO,), the boric acid (H3BO3) ionization, and the shear viscosity of seawater [10, 
11]. Sound absorption is quantitatively described by the absorption coefficient a with 
units of dB/m. One of the most accurate expressions for the frequency dependent absorp- 
tion coefficient is provided by Fisher and Simmons [10] and is based on experimental 


data; 





a. [dB/m] = es re au ES fos (2.7) 
Tite deeT 
The terms A, B, C, depend on temperature and hydrostatic pressure whereas the relaxa- 


tion frequencies f, and f, are functions of temperature only. Specifically, after the ap- 


proximation given in [10] 


f, [Hz] =780e"”, (2.8) 
fy [Hz] = 42000e""", (2.9) 
A=8.3x10° (§/35) rere er (2.10) 
B=0.022(S/35)e7"'*7'%, (2.11) 
and 
CHARM eras (2.12) 


where Z is the water depth in km, T is the temperature in °C, S is the salinity in ppt and 
pH is the acidity. The absorption coefficient in seawater and frequencies between 100 Hz 
and 1 MHz is shown in Figure 2 (from Equation (2.7) with parameter values given in the 
figure caption). In shallow water (Z < 150 m), depth, temperature and salinity have little 


effect on the absorption of sound. 





o. (dB/km) 


I bo Pog | | I ) Paaii 
ISISISIA BS RP RAP eS BAHASA PAR eee ae as 


= SSIS AO Se ee Pea ae oe 
1=1=1 
I 1 1 














Frequency (Hz) 


Figure 2. Absorption Coefficient (dB/km) in Seawater (pH = 8, S = 35 ppt, Z= 0.1 
km, T= 15°C). 


Now that simple forms for energy spreading and sound absorption are defined, a 


general approximation of the transmission loss in shallow water is given by [9, 11] 


20logr +ar rar 
TL [4B ref. 1m] = TL gpeading FTL absorption = 10logr+10logr,+ar r>r,’ ee) 


where the transmission range r and the transition range 7, are in meters and the absorp- 


tion coefficient @ in dB/m. The range and frequency dependency of TL for a fast bottom 
environment is illustrated in Figure 3 (for the parameters given in the captions of Figures 
2 and 3). The frequencies chosen to evaluate Equation (2.13) are the lowest and highest 
operating frequencies of Seaweb, f; = 9 kHz and f> = 14 kHz, respectively. A third fre- 
quency (f3; = 50 kHz) much higher than /> is used to demonstrate the logarithmic increase 
of the TL at high frequencies. In Figure 4, the TL is plotted as a function of frequency for 


the maximum transmission range (5 km) according to Seaweb specifications. 
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Figure 4. 


As seen in Figure 4, for a frequency of 12 kHz (the mean frequency of Seaweb), 
the TL is approximately 45 dB; whereas for a frequency of 50 kHz the TL increases loga- 
rithmically assuming the value of 110 dB. The increase of the TL at high frequencies is 
the reason the available communication bandwidth of Seaweb is limited to 5 kHz. Un- 
derwater communication channels are severely band limited compared to the channels 
used in RF communications. Because of the limited bandwidth, it is very difficult to 
achieve a reliable link at high bit rates. By increasing the data rate, the period of the mes- 
sage bit decreases and the bandwidth of the signal increases. For the sake of comparison, 
each channel of the 802.11g wireless network according to the IEEE standard has a 
bandwidth of 25 MHz and a maximum throughput of 54 Mbps, whereas an underwater 
network can only have a few kHz of bandwidth and a much lower data rate. 

B. AMBIENT NOISE 

The underwater environment everywhere contains ambient noise. Ambient noise 
is the total noise background minus the noise generated by the transmitting equipment 
used in a certain application. Noise observed in the ocean exhibits strong frequency de- 
pendence as well as geographical dependence. In deep water, there are many sources of 
ambient noise and each source dominates in a different frequency band. At low frequen- 
cies (below 20 Hz) the main noise sources are the ocean turbulence and the various tec- 
tonic processes such as earthquakes and volcanic eruptions [13]. In the frequency band 50 
to 500 Hz, distant ship traffic and distant storms are the main sources of noise [12]. Be- 
tween 0.5 and 50 kHz, the detected noise is associated with the breaking of the sea sur- 
face and bubble formations. The state of the local sea surface is weather dependent and 
closely connected with the prevailing wind speed. The greater the wind speed the greater 
the ambient noise contribution. Above 50 kHz, thermal noise due to the molecular agita- 
tion of the water molecules predominates [10]. The interested reader can read more about 
ambient noise and its sources in References [14] and [15]. The ambient noise spectrum 
level (ANL) due to the combined effects of the aforementioned sources is illustrated in 
Figure 5. In the frequency band (9 to 14 kHz), at which Seaweb operates, the wind speed 
(or sea state) is the main source of noise. For example, at 12 kHz and a wind speed of 15 


knots, the ambient noise spectrum level is approximately 45 dB ref. 1 Pa. In general, the 
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ANL decreases when the frequency increases, and falls at a rate of 15 dB/decade when 


the frequency range is between 0.5 and 50 kHz. 
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Figure 5. | Average Deep-Water Ambient Noise Spectrum Level [From Ref. 12]. 


Wenz curves [15] (or modified ones like Figure 5) are commonly used to estimate 
the ambient noise spectrum level under various conditions (e.g., shipping density, sea 


state, etc.). The following expressions [16] are a quantitative description of Wenz curves: 


ANL parbutence = 17 — 30 log(f), (2.14) 
ANL ghipnine = 40 + 20(D — 0.5) + 26 log( f) — 60log( f + 0.03), (2.15) 
ANL wing = 50+ 7.5(w"* )+ 20 log(f)— 40 log(f + 0.4), (2.16) 
ANL thermal = ~15+ 20log(f), (2.17) 


and 


ANLjurbutence ANL Shipping ANLyind ANL therm 


ANE an “FOL 1 © 410 © 410 © 410 =| (2.18) 








where D is the shipping density with values between 0 and 1 (1 for heavy traffic), w is the 


wind speed in m/s, f is the frequency in kHz, and ANL is the ambient noise level 


Deep Water 
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in deep water with units of dB ref. 1 uPa. Equation (2.18) is evaluated for various ship- 


ping densities in Figure 6 and for various wind speeds in Figure 7. 
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Figure 6. 
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Ambient Noise Level in Deep-water for Wind Speed w = 15 m/s. 
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In shallow water, ambient noise levels are very difficult to define because two ad- 
ditional noise sources, man-made noise and biological noise (dolphins, snapping shrimp, 
etc.) exist. These sources vary greatly both in time and location. Man-made noise is 
higher close to harbors and coastal areas with industrial activity. Noise from biologics is 
in general unpredictable. It is quite difficult to estimate beforehand and existing meas- 
urements (if any) are only valid in the specific location they were taken. However, in the 
absence of industrial and biological noise, the ambient noise level depends mostly on the 
strength of locally generated wind and is about 5 dB higher than the deep water levels for 
the same frequency [15]. 


Both transmission loss and the ambient ocean noise contribute to the degradation 
of the acoustic signal as it propagates away from the source. By combining Equations 
(2.13) and (2.18) and assuming that the noise generated by Seaweb equipment is zero, we 
can calculate the signal-to-noise ratio (SNR) of the underwater communication channel 


from 
SNR [dB] = SL+DI, +DI, -TL—ANL, (2.19) 
where SL is the source level in dB ref. 1 Pa, 1m, and DI,, DI, are the directivity indi- 


ces of the transmitter and the receiver, respectively. The directivity index is a quantitative 
measure of the focusing of acoustic energy, and for an omni-directional source or re- 
ceiver, DI=0. The SLof the low frequency omni-directional “telesonar” transducer AT- 
A08 (a candidate for Seaweb) according to the OEM specifications [17] is 180 dB ref. 1 
uPa, 1m at full power. The resulting SNR of the channel is illustrated in Figure 8. We can 
see from this figure that for a frequency f= 14 kHz and a transmission range of r = 5 km, 
the SNR is approximately 76 dB. Needless to say, this is only a rough and rather optimis- 
tic estimate of the SNR found in a real underwater environment. The frequency depend- 
ent SNR for various transmission ranges is shown in Figure 9. The reason the frequency 
band of 9 to 14 kHz is chosen for the Seaweb network is clear since this band accounts 
for the flat portion (SNR is almost constant) of the curves corresponding to r = 3 and r = 


5 km. 


13 


0.5). 





16000 











14000 











8000 10000 12000 
Frequency (Hz) 


6000 


4000 





2000 


Q 
° a 
w |! r= 
NZ /I 
mx], Va) 
~~) imal 
HO WH---- o 
1 3 
| o 
Q 
n 
cl pa as Ragga ell grenh Ms Syeleeen, e l tan d  I oe,o8 = 
= 
— 
Es 
— 
are n 
o 
“— 
€ 4a 
ee ae ee, eee ere eee eee g 38 
Q ion 
fo?) 
S ia 
o Coy 
Pee ae Se Se a Se See Se Se Sit SEAS ae a ES a eS ee Cc lo?) 
2 
g 
= ise] 
2 o 
© Q 
F © 
OQ 
o 
3 
3 
o 
N 
o 
s 
u 
& 
~ 
Z 
N 
































Figure 8. 


0.5). 


SNR for Different Transmission Ranges (wind speed 15 m/s, D 
14 


Figure 9. 


C; TEMPORAL AND SPATIAL VARIABILITY OF THE CHANNEL 


The characteristics of the underwater communication channel vary both in time 
and location. The transmitted signal is subject to multipath propagation and/or Doppler 
spreading. Together, these factors may cause a severe degradation in the received signal 
strength. 

1. Multipath Propagation 

The acoustic wave as it travels through the ocean rarely follows a single path. Due 
to reflection at the sea boundaries, refraction, and volume scattering within the water col- 
umn or the sub-bottom structure caused by random inhomogeneities, multiple echoes of 
the transmitted signal arrive at the receiver. These time-delayed echoes vary both in am- 
plitude and phase and interfere with one another. This effect is called multipath propaga- 
tion or time spreading since travel time is different for each path. Multipath characteris- 
tics are different in shallow and deep-water channels. In shallow water, the multipath is 
severe because of repeated reflections from the sea surface and bottom. Each reflection 
causes an attenuation of the sound wave prohibiting long-range communications without 
a significant loss of acoustic energy [13]. The relative delay time between the first and 


i-th arrival of the transmitted signal at the receiver is called excess delayz,. The maxi- 


mum excess delay 7,,,. between the first arrival and the arrival of the last distorted rep- 


lica of the original signal is much greater in an acoustic communication channel than in a 
land-based RF channel. This is mainly because sound in water travels slower (~ 1500 
m/s) than the electromagnetic waves travel in the air (3x10* m/s). In the frequency do- 
main, the bandwidth over which the channel passes all spectral components with ap- 


proximately equal gain is called the coherence bandwidth B.[18] and is a useful tool for 


comparing different multipath channels: 


1 
Bo = —. (2.20) 
2. Doppler Spreading 
The underwater environment is time variant. The medium itself is constantly 
moving due to surface and/or internal waves and the source and receiver platform may 


also move. In shallow water, surface scattering is anticipated to be the most important 
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contributor to the overall time variability of the channel. Surface scattering is caused by 
the weather dependent roughness of the sea surface. The frequency of an acoustic wave 
incident upon the moving surface of the ocean is affected by spectral broadening. This 
phenomenon is called Doppler spread. The Doppler spread for the case of a sound wave 


of frequency / reflected at the sea surface is given by [13] 


3/2 
B, -(u.7ssa0?)[ 2282) (2.21) 


Cc 


where B, is the Doppler spread in Hz, w is the wind speed in m/s, @ is the grazing 
angle, and c is the speed of sound ( ~ 1500 m/s). For a transmission range much greater 
than the channel depth (z >>7 ), @ is very small andcos@ ~1. In Figure 10, the Doppler 
spread for a frequency f=14kHz and wind speed w=15 m/s (sea state 3) 
is approximately 9 Hz and the frequency spectrum of the reflected wave 


is f +B, =14+0.009 kHz. 
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Figure 10. Doppler Spread Due to Surface Reflection. 
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The second parameter used to describe the time variability of the channel is the 
coherence time 7/., which is inversely proportional to the Doppler spread (for simplicity 


the Doppler spread is considered to be the only factor contributing to coherence time), 
T.=—. (2.22) 


The meaning of coherence time is that two signals arriving at the receiver with a time 


separation less than 7,, will be affected by the channel in the same way [18]. 


The relation between the aforementioned parameters (,,,., B., Bp, T) that de- 


scribe the overall variability of the channel leads to four distinct types of fading that are 


summarized in Figure 11: 






Fading Based on Coherence Bandwidth 


Frequency Non-selective Frequency Selective Fading 


Fading Based on Coherence Time 
Fast Fading Slow Fading 


Figure 11. Different Types of Fading Channels [After Ref. 18]. 









The first two types of fading are frequency selective and frequency non-selective 
fading. A frequency non-selective channel implies that the coherence bandwidth of the 


channel B. is much greater than the bandwidth of the transmitted signal B, or equiva- 


lently the maximum excess delay r,,,. is much less than the digital symbol (bit) pe- 


IX 


riod 7, . In this case, which is the best for wireless communications, the channel induces 


no intersymbol interference (ISI). On the other hand, when frequency selective fading 


(B. < B,) occurs, then 7, < r,,,, and ISI is unavoidable. 


The other two types of fading are based on the coherence time of the channel. The 
first one is called a “‘fast-fading” channel and the second one a “slow-fading” channel. In 
17 


a fast-fading communication link the bit duration is greater than the coherence time of the 
channel. The channel impulse response changes rapidly during the time it takes to trans- 
mit one symbol. In this type of fading a severe distortion of the signal occurs and a reli- 


able communication link cannot be maintained. In a slow-fading channel, 7, << Ti, 
(B, << B,), and less distortion is introduced as a result of the wave propagation in the 


water. 


In this chapter, we have seen the reason the shallow-water channel is considered 
to be one of the most challenging channels for wireless communications. The system en- 
gineer must take into account all the limiting factors (TL, ambient noise, fading) imposed 
by the channel in order to design a reliable communication link. Because of the time- 
varying multipath, there is always a trade-off in choosing the transmission rate. A higher 
bit rate implies greater bandwidth, higher transmission loss and ISI but an almost negligi- 
ble channel variation over a single-bit interval. In order to undo the effects caused by the 
channel, a statistical model of the channel impulse response is needed. In the next chap- 
ter, we describe the steps needed to extract the channel impulse response by using one of 
the existing channel models called MMPE, which is based on the parabolic approxima- 


tion of the wave equation. 
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HI. MODELLING OF THE UNDERWATER CHANNEL 


In order to allow for the implementation of an efficient communication scheme, a 
statistical representation of the underwater channel is needed. In underwater acoustics 
various numerical methods are used to describe the propagation of sound. All current 
methods solve the wave equation for a specific environment by following two different 
theoretical approaches. The first approach is based on wave theory (normal-mode, para- 
bolic equation models) and the second one on ray acoustics (ray models) [12]. Each the- 
ory has its own advantages and disadvantages. Wave theory models give a complete solu- 
tion and are best suited for shallow-water and low-propagating frequencies. Ray models 
provide a practical visualization of the propagation of sound in space and give their best 
results at high frequencies [11]. 

A. MONTEREY-MIAMI PARABOLIC EQUATION (MMPE) MODEL 

In this thesis the channel impulse response is predicted using the Monterey-Miami 
parabolic equation (MMPE) model [7]. The MMPE is a range-dependent, wave- theory 
model based upon the parabolic approximation of the wave equation. Since the purpose 
of this research is not the model itself, only an outline of its theoretical foundation is pre- 
sented. We begin by representing the time harmonic acoustic pressure field in cylindrical 


coordinates assuming azimuthal symmetry, 
P(r,z, ot) = p(r,zje. (3.1) 


Next we substitute Equation (3.1) into the linear, homogeneous wave equation in cylin- 


drical coordinates, 


2 


oe z)=0, (3.2) 


V’ p(r,z)— 
P(r) c(r,z) Ot? 





where c [m/s] is the speed of sound in water and V* the Laplacian defined as 


Vo (3.3) 
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This leads to the time-independent Helmholtz equation, 


2 
@O 


c(r,zy 


which can be factored by introducing the operator Q,, defined as 


0, =VH+eEFt1, (3.5) 





V’ p(r,z)+ P(r,z) =9, (3.4) 


with 
e=n -l, (3.6) 
n=, (3.7) 
Cc 
and 
1 0 
=, 3.8 
Mt eo (3.8) 


In the above expressions, is the index of refraction, c, is the reference sound 
speed of the ocean volume (usually taken as 1500 m/s), and k, =q@/c, is the reference 
wavenumber. Note that the coordinate dependency of n and c is dropped for simplicity 


but is always implied. 


If we take into account the cylindrical spreading, which is typical in shallow water 
and properly factor the Helmholtz equation, the out-going acoustic pressure may then be 


defined as [19], 


R = ikyr 
P(r,2) = Py pw (r.ze, (3.9) 


where P, is the pressure amplitude at distance r= R,, and y(r,z) is the envelope func- 
tion, or PE field function. After substituting Equation (3.9) into (3.4) the parabolic equa- 


tion for the PE field function is formulated: 


fe) : . : 
oa ~ “ky + ik,Q,,W = —tky HW, G3. 10) 
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where 


H, =1-Q, (3.11) 
is an operator, defining the evolution of the PE field function in range. By integrating 


Equation (3.10) and assuming that, is constant over a distance Ar, the relationship 


between the values of y at different ranges may be defined as 


w(rt+Ar) =e nr Oy(r), (3.12) 


pion aac propagator that accounts for the change in the value of the PE field 


where e 
function y with range. This propagator is computed using a split-step Fourier (PE/SSF) 
method described in [20]. The discrete fast Fourier transform (FFT) subroutine allows the 


PE/SSF implementation to be represented by [7] 


~iky 


y(r+Ar,z)=e 


Ar 
2 


r Ar 
op (r+ Ar,z) -ikyArT.,, (k, thy Uy (12) 
ner FrT|e tua are mews vie} (3.13) 


The U,, and Ts operators are defined as 


U,, =-(n-]), (3.14) 


and 


2 
- k 
T, =1-,/1-| =], 3:15 
; [F (3.15) 


where k,is the vertical wavenumber. The output of the model is the PE field function y 


from which other acoustical quantities can be computed. For example by substituting 


Equation (3.13) into (3.9) we obtain the ocean frequency response, 


G(f) = p(r,2) =(rz)e™ (3.16) 


For broadband signals the MMPE model must be run for each frequency in the 
band in order to compute the frequency response at a given range and depth. In the fol- 
lowing sections of this chapter, the theory behind the various perturbations modeled in 


this thesis is described and the results of the simulations are presented. 
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B. INTERFACE ROUGHNESS PERTURBATION 
In order to define the interface roughness, we assume a 2-D interface spectrum of 


the form [21] 


Wo (k,)=——,, (3.17) 
(142. k?)? 


corr r 


where k, is the horizontal spatial wavenumber vector given by 


k,=VK?4+P, (3.18) 


K and L are the horizontal wavenumbers in the x- and y-directions, respectively, Leo 1S a 
correlation length scale, f is the spectral exponent and yz is a normalization factor de- 


fined in terms of the root-mean-square (rms) roughness o [22], 
u -2(£-1Jo%2.,. (3.19) 


For the determination of the interface roughness affecting forward propagation, 


the one-dimensional (1-D) spectrum along the x-axis W;” (K ) is required. This is com- 


puted by taking the 1-D transform of the 2-D interface spectrum WL(K,) =", (K 2b) 


along a slice at y =0, defined by 
Wi? (K)=[_ WS (K,L)dL. (3.20) 


In cylindrical coordinates Equation (3.20) yields 


Bil 
(+2 K?) ? 2, (3.21) 


corr corr 


lf 


and I is the Euler gamma function. In order to illustrate the interface roughness pertur- 


W(K) = yo L 


where y is given by 





(3.22) 


bation to the underwater environment, the TL vs. depth and range is computed for the 
frequency f= 11.5 kHz, a bottom depth of 150 m, and a source depth of 75 m. The results 


are shown in Figure 12 for an rms roughness value of 2 m. By looking at the two plots 
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two observations can be made. For the case of the unperturbed environment (top plot) the 
down-refracting sound speed profile causes significant interactions with the seafloor. The 
upper figure also displays a clear interference pattern due to the coherent summation of 
modes. On the other hand, when the interface roughness is added, the bottom interactions 
with the rough boundary create a diffuse field due to the incoherent scatter of energy and 


the breakdown of coherence in the mode interference pattern. 
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Figure 12. Interface Roughness Perturbation Compared to the Unperturbed Environ- 
ment (Source Depth = 75 m, f= 11.5 kHz, Range 0 to 1 km). 
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C. BOTTOM VOLUME PERTURBATION 

1. Volume Sound Speed Fluctuations 

For the computation of the volume sound speed fluctuations affecting the forward 
propagation, the two-dimensional (2-D) vertical spectrum of the volume fluctuations 
wy ky) is required, where k, =/K°+M° is the wavenumber in the (x,z) plane. This 
can be determined by first modeling the sediment volume sound speed perturbations with 


help of the three-dimensional (3-D) volume spectrum WYK ,L, M) defined by [23] 


Wi(K,L,M) ~ PAB (K: +B)smey a, (3.23) 
TW 


where M is the vertical wavenumber, # is the spectral strength constant, and A is the 
horizontal-to-vertical aspect ratio used to describe the anisotropy of fluctuations in the 


sediment. 


The 2-D vertical volume spectrum wy ) (K,M) in the (x,z) plane is computed by 
integrating (3.23) over Las follows [22]: 
B4+2 


W!))(K,M)= f W!)(K,L,M)dL = BABI (K? sees) dL. (3.24) 
T 


= 0 


For the specific case of @ =2, A~5, and Bx5x10~ [23], Equation (3.24) yields [22] 


N| vw 


WS) (K,M) = | W{O(K,L,M) dL =a'[ 25K? +M?] ?, (3.25) 


—00 


where a’ =1.25x10>. 
2. Volume Density Fluctuations 


The variability in bottom volume density p is treated through the effective index 


of refraction n' defined by [7] 


2 
nn? =n + a ly s a ; (3.26) 
2k, | P 2\ p 


where k, is the reference wavenumber and n is the index of refraction. By considering 





the fact that the sediment properties are horizontally stratified, and assuming that the en- 
vironment is range-independent over a range step Ar , Equation (3.26) simplifies to 
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2 2 
per eee Lay Bf 2 (3.27) 
2k,| p Oz 2\p o& 


The definition of the index of refraction in the bottom is 


i= i. (3.28) 
C, 
where 
Cy = Cy, (1+b,+6,)= C, + 0C,. (3.29) 
The factor dc, 1s given by 
OC, =C, (b,+6,), (3.30) 


where c, is the mean bottom sound speed at the interface, 6, is the zero-mean random 


perturbation, and b, = g / c,, 1s the normalized gradient of the bottom sound speed. 


The relative fluctuations in bottom volume density are related with the relative 


fluctuations in sound speed through the expression [23] 





Sp _ “{e.=Pal 22) (ay) % (3.31) 
Py 22) - P, \ a 
where 
Cy = Cy, (1+8, ), (3.32) 
and 
yp oe Bove, (3.33) 
2p — Pp, 


In this case, p, = 2650kg / m° is the density of the grain, and ,, c, are the average val- 


ues of the density and sound speed in the sediment, respectively. From Equation (3.31), 


we obtain for the bottom-volume density with perturbation, 
Oc Oc 
P=p,top= Po* Po(27)— = Pp 1+27—|. (3.34) 
0 0 


The first and the second partial derivatives of (3.34) with respect to depth 
z (depth gradients in p, or c, are neglected) yield 
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2 FUE ia, (3.35) 
OZ Gy 


and 





2 2 
8 221s O52), (3.36) 


Gz? C, 2° 
These are then incorporated into the definition of the effective index of refraction, Equa- 
tion (3.26), within the bottom volume. The bottom-volume perturbation in terms of the 
TL vs. depth and range is displayed in Figure 13. As shown in this figure, the scattered 
acoustic energy at short distances from the source is penetrating into the bottom. How- 
ever, this has only a minor impact in the resulting interference pattern compared to the 


one created by the unperturbed environment. 
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Figure 13. Bottom Volume Perturbation: TL vs. Range and Depth (Source Depth = 
75m, f=11.5 kHz, Range 0 to 5 km ). 
D. TURBULENCE PERTURBATION 
Another mechanism affecting the coherence structure in an underwater environ- 
ment is the small-scale fluctuations of sound speed within the water volume caused by 
turbulence. The turbulent perturbation field is approximated as a random realization of 


perturbations that are based on an accepted power spectral density of the turbulence. 
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The variance of the perturbations to the index of refraction (7) is given by [22] 








var( (7) = J «_8,(K,L,M)dK dLdM = ee i S(kK)d*k, (3.37) 


1 
(27) 
where S', (k) =S'(K,L,M) is the spectral density of the perturbations and k= (K ,L,M ) 


is the 3-D wavenumber vector. If the variance is normalized according to 





var(u(F)) = (lu) ) = [we @a’r, (3.38) 


me 


where ( ) denotes averaging , then by combining Equations (3.37) and (3.38) we obtain 


aul 





I, Uw (F)d*r = of S,(k)d*k. (3.39) 


Now, if F, (k) is the perturbation transform, the index of refraction perturbation (7) can 


be written as 


=) 1 f Ty ikeF 73 
HA = Ts J F (ke""d’k, (3.40) 


and its spatial autocorrelation is given by 





j UP) (F -7)d?r = oni {|} i F (hed «[j [ Baer a? ila (3.41) 


For r’=0, Equations (3.39) and (3.41) yield the desired relation between the signal spec- 


trum F (k) and the 3-D spectral density of the perturbations S$ Ak) [22]: 
F (| =| volumex S,,(%)] 3.42 
| me )| =[ volume x mi )| , (3.42) 
Similarly, along the 2-D plane of propagation (F = (x,0, z)) , the relation between the 2-D 
spectral density V,(K,M), and the 2-D spectrum of the perturbation field G,(K,M), 
can be written as [22] 
1/2 
|G,(K,M)|=[AreaxV,(K,M)] (3.43) 


The turbulent perturbation field is approximated through the use of the method of 
structure functions [24], in which large-scale perturbations are not accounted for. The 


three-dimensional spectral density function E(k), can be expressed as [25] 
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E(k) = Ae??k""?, (3.44) 
where k7''* = (K 24D +M? Vie , € is the energy dissipation rate, and A is a scalar mul- 


tiplier. In order to prevent E(k) from growing unbounded when the value of k is small, 


a lower wavenumber threshold k, is defined so that the k~'"’’ dependence is [22] 


1 


Also, by applying the high-pass filter 
k? 


the 3-D spectral density function E(k) is forced to assume the value of zero when k =0. 


Subsequently, a high-frequency cutoff is established expressed as [25] 


R,(k) = exp (4) } (3.47) 


where R,(k) is the Batchelor spectrum, g is an order unity factor, and k, is the Batchelor 


bs =[ =) , (3.48) 


2 
V Kp 


wavenumber defined by [22] 





In the above expression, € [W/kg] is the depth averaged kinetic energy dissipation rate, 


x, ~ 10°’ m’/s is the thermal diffusivity of the sea water, and v=1.4x10~% m’/s is the 


kinematic viscosity. In this theoretical development, € is given by [26] 


r=%3([S amet ans cosh"! (#2) (3.49) 
ff, \La Si ¥ 


where f, [rad/s] is the inertial frequency, bDE,,, ~ 0.5m is a measure of internal wave 
intensity, j, is the mode number, and N(z) [rad/s]is the buoyancy frequency expressed 


in terms of density p(z) and potential density p,(z), 


ae 1 ap, a 1/2 
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By following the theoretical approach described above, the three-dimensional spectral 
density S'(4;z) can now be expressed by [22] 
A, k? 


-2/3 


where 4, [m~*] is the turbulent strength parameter. In this implementation, 4, will be 


adjusted manually to produce a specific desired rms fluctuation. 


The desired 2-D spectral density is now formulated by integrating over L (y =0) 


only the portion of the 3-D spectral density that varies with the turbulent spectrum. This 
portion is given by [22] 














1 
Thus, the 2-D spectral density ®, ,(k) can be derived from 
1 f 1 1 ¢ 1 
®,,(K,M)= an J (ke +e dL = on J (K?+2+M +key dL 
1 ¢ 1 1 *F 1 
= on | (2+Ey"s dL = jee | py aL, (3.53) 
. ; [4 | 
a 
where a” = K* +M’* +k°*. With help of integral tables, Equation (3.53) yields 
1 \vzr@) 
© p(Kt)=| 5 re (3.54) 


By reapplying the terms that do not vary with the turbulent spectrum, the resulting 2-D 
spectral density ®,,,(k) is defined by [22] 


rtbiey=|[ S| Ay = Jv (3.55) 


2n NIT (2) (C4 Le +e 


The perturbation to the underwater environment because of turbulence is illustrated in 





Figure 14 for an rms value of 2 m/s. The breakdown in the coherent structure of the 
modes is obvious, resulting in a highly diffuse field with random patterns of ray-like 


propagation (yellow, high-energy paths in Figure 14). 
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Figure 14. Turbulence Perturbation: TL vs. Range and Depth (Source Depth = 75 m, 
f= 11.5 kHz, Range 0 to 5 km). 


E. INTERNAL WAVE PERTURBATION 

Larger-scale fluctuations in sound speed may also be generated from internal 
waves. The internal wave perturbation accounts for the vertical movement of water 
masses inside the sea body due to buoyancy forces. In this thesis a combination of a lin- 
ear (sinusoidal) and a nonlinear soliton internal wave perturbation to the sound speed is 


modeled. The linear internal wave perturbation is expressed as a sum of sinusoids [27], 


5 
Ac, (7,2) =C exp(-2) >) cos(K,r), (3.56) 


i=l 
where Ac (7.2) is the sound speed perturbation at a given range r and depth z, 
B=25m, K,=2z/4,, 2, = 1500-300(i-1) m, and C is defined in such a way that 


the maximum value of Ac,,, (42) is 1.5 m/s. 


The nonlinear soliton internal wave perturbation Ac,,, (r,z) is defined as [27] 


2 
6 — 
Aes (raz)=CZexp(-2 94 o( 8S "| ‘ (3.57) 





where B = 25m, A, =10exp(-0.3(i-1)), R, = 4800m, R,,, = R,,-200(7-i) m, 
D, = ,{34300/ A, m, and C is defined such that the maximum value of the sound speed 
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perturbation is 12.5 m/s. The internal wave perturbation in terms of the TL vs. depth and 
range is shown in Figure 15 (range 0 to 1 km) and Figure 16 (range 4 to 5 km). At short 
ranges the perturbation is mostly due to the linear internal waves and the modes seem to 
remain coherent. At long ranges from the source (Figure 16) the nonlinear soliton internal 
waves cause a much stronger perturbation to the field resulting in a more random like 


dispersion of the energy. 
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Figure 15. Internal Wave Perturbation Compared to the Unperturbed Environment 
(Source Depth = 75 m, f= 11.5 kHz, Range 0 to | km). 
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Figure 16. Internal Wave Perturbation Compared to the Unperturbed Environment 
(Source Depth = 75 m, f= 11.5 kHz, Range 4 to 5 km). 


F. DOPPLER PERTURBATION 


The influence of Doppler due to source/receiver motion is of great interest in un- 


derwater acoustic communications. Because of the Doppler effect the frequency of a 


transmitted signal suffers from spectral broadening. The actual frequency /, (0) , propa- 


gated from a source that is moving with speed v,, can be expressed as [28] 





f.(0)= Sr ah f1+[2)oos(o-#,)} (3.58) 
1-( 2 Jeos(#-4,) . 
Cc 
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where /, is the transmitted frequency, c is the sound speed, @ is the angle of propaga- 
tion (relative to the horizontal), and ¢, is the angle of source motion. As the sound trav- 


els through the medium it follows multiple paths. Therefore, even a single transmitted 


frequency will result in multiple received frequencies at the receiver. 
If the bandwidth BW, of the transmitted signal is 
BW, = frmax — Fria » (3.59) 
where fr min 20d fp max are the minimum and maximum frequency components of the 


transmission, then the Doppler shifted bandwidth for forward propagation is given by 


[28] 


BY = BW, +©( fia foie 
Cc 





sin(¢,)|), (3.60) 





or 


BWs- = BW, +E (Firs 
Cc 





sin(¢s)|+ Srmin)- (3.61) 


Equation (3.60) accounts for a source motion toward the receiver (7 + ) whereas (3.61) is 


for a source motion away from the receiver (7 -—). 


If the receiver is moving also, a similar theoretical approach yields [28] 








fe D 
BW* =BW, | F osmax + fonin (Sin ($e )|), (3.62) 
or for the case when the receiver is moving toward the source 
if D 
BWi =BW, | Frmax |Sin ($e )|+ Fomin ) (3.63) 





where v, is the receiver speed, ¢, is the angle of receiver motion, fp nin ANd fo max ALE 
the minimum and maximum frequency of propagation, and BW, is the bandwidth of 
propagation given by 

BW = Spies J Paint (3.64) 
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G. RESULTS 

The results of the various perturbations described above are presented in this sec- 
tion. The specific parameters defining the characteristics of the underwater environment, 
the source/receiver, and the resolution of the calculations are summarized in Table 1. For 
each simulation the MMPE model was run over a 2048-Hz bandwidth (1024 frequen- 
cies), and the channel impulse response was calculated for a source-receiver distance of 5 
km and a source depth of 75 m as in a typical Seaweb setup. For all modeled perturba- 
tions, the transmission loss is presented as a function of depth and frequency (Figure 17 
through Figure 23 - top plot). The magnitude of the channel impulse response at a depth 
of 75 m versus the travel time of the acoustic wave is illustrated in Figure 17 through 
Figure 23 (bottom plot). The impulse responses are normalized by the magnitude of the 


unperturbed response such that its values lie between —land 1. 


From the results presented in the following figures some interesting observations 


can be made. The typical excess delay r,,,, between the first and last arrival of the trans- 


mitted signal is about 0.15 seconds (150 ms). For a symbol duration T. =0.42ms _ corre- 


sponding to a chipping rate of 2400 chips-per-second (used in the DS-SS implementation 


described in Chapter V) channel-induced ISI will occur (T. < r,,,,. ). 


Among the individual perturbations examined, the one corresponding to the inter- 
face roughness (Figure 18) seems to be generating the worst fading environment. This 
was expected since, as already mentioned, in shallow-water channels the interface rough- 


ness is typically the dominant scattering factor. 
Figure 23 represents the worst-case scenario where all perturbations are combined 
and Doppler (vs = 15knots ~ 8m/ s) is added to account for the spectral broadening in- 


duced by the channel. The effect of Doppler is obvious by comparing the results illus- 
trated in Figure 22 (all perturbations and no Doppler) and Figure 23. It is obvious that the 
underwater channel is far more challenging for communications when Doppler is in- 


cluded. 


34 


















































Parameter Value Remarks 
Main Parameters 
Number of depth points saved 3277 
Minimum depth saved 0m 
Maximum depth saved 200 m 
Number of range steps saved 10 
Minimum range saved 0km 
Maximum range saved 5 km 
Computational range step size 0.4m 
Maximum computational depth 500 m 
Number of computational depth 16,384 Radix-2 integer required for FFT 
points 
Reference sound speed 1500 m/s 





Source Data 











Source depth 75m 

Center frequency 11.5 kHz 

Frequency bandwidth 2048 Hz 

Number of Frequencies 1024 Radix-2 integer required for FFT 





Sound Speed Profile Data 





Water column sound speed 


1500 m/s 


(unperturbed) 
Bathymetry of the Water-Bottom Interface 


Isospeed 





Bottom depth (unperturbed) 


150m 


Flat 





Acoustic Parameters of the 


Medium below the Water-Bottom Interface 

















Bottom sound speed 1600 m/s 
Sound speed gradient I/s 

Relative density (bottom/water) 1.6 
Compressional attenuation 0.1 dB/km/Hz 
Shear speed 0 m/s 

Shear attenuation 0 dB/m/kHz 





Acoustic Properties of the Deep 


Layer 











Sub-bottom sound speed 2000 m/s 
Sound speed gradient I/s 
Relative density (sub-bottom/water) 3 g/em* 





Compressional attenuation 


0.25 dB/km/Hz 





Shear speed 


0 m/s 








Shear attenuation 





0 dB/m/kHz 
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Parameter Value Remarks 
Perturbation Parameters 
Source speed 15 knots (8 m/s) Doppler perturbation 





Bottom interface roughness (rms) 


2m 


Bottom interface perturbation 





Bottom volume sound speed vari- 
ability (rms) 


15 m/s 


Bottom volume perturbation 
(with corresponding density per- 
turbation, as described in text) 





Internal wave (linear) sound 
speed variability (length scales & 
max amplitudes) 


Internal wave (soliton) sound 


1500 m & 1.5 m/s 
1200 m & 1.5 m/s 
900 m & 1.5 m/s 
600 m & 1.5 m/s 
300 m & 1.5 m/s 


4800 m & 12.5 m/s 


First part of IW perturbation (al- 
ways combined with second part 
of IW) 


Second part of IW perturbation 


speed variability (ranges & max 
amplitudes) 


3800 m & 9.3 m/s 
3000 m & 6.9 m/s 
2400 m & 5.1 m/s 
2000 m & 3.8 m/s 
1800 m & 2.8 m/s 


(always combined with first part 
of IW) 

















Water column turbulence (length | 100m & 2 m/s Turbulence perturbation 
scale and rms) sound speed vari- 
ability 

Table 1. .MMPE Parameters. 
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Figure 17. No Perturbation: TL as a Function of Depth and Frequency (top plot) and 
Normalized Impulse Response at Depth 75 m vs. Travel Time (bottom plot). 
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Figure 18. Interface Roughness Perturbation: TL as a Function of Depth and Fre- 
quency (top plot) and Normalized Impulse Response at Depth 75 m vs. Travel Time (bot- 
tom plot). 
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Transmission Loss (dB re 1m) - Bottom Volume Perturbation 
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Figure 19. Bottom Volume Perturbation: TL as a Function of Depth and Frequency 
(top plot) and Normalized Impulse Response at Depth 75 m vs. Travel Time (bottom 


plot). 
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Transmission Loss (dB re 1m) - Internal Wave Perturbation 
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Figure 20. Internal Wave Perturbation: TL as a Function of Depth and Frequency 
(top plot) and Normalized Impulse Response at Depth 75 m vs. Travel Time (bottom 


plot). 
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Transmission Loss ( dB re 1m) - Turbulence Perturbation 
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Figure 21. Turbulence Perturbation: TL as a Function of Depth and Frequency (top 
plot) and Normalized Impulse Response at Depth 75 m vs. Travel Time (bottom plot). 
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Figure 22. All Perturbations Combined (no Doppler): TL as a Function of Depth and 
Frequency (top plot) and Normalized Impulse Response at Depth 75 m vs. Travel Time 
(bottom plot). 
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Transmission Loss (dB re 1m) - All Perturbations (with Doppler) 
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Figure 23. All Perturbations Combined (with Doppler): TL as a Function of Depth 
and Frequency (top plot), and Normalized Impulse Response at Depth 75 m vs Travel 
Time (bottom plot). 


This chapter presented a brief description of the MMPE model and an outline of 
the theory behind the various perturbations modeled. The MMPE model was used to pre- 
dict the channel impulse response for each individual perturbation so that the implemen- 
tation of the DS-SS scheme with the CRV updating algorithm could be evaluated in 
Chapter V. In the following chapter the CRV decomposition and its respective updating 
algorithm is presented, and the reason that the CRV can be used as an alternative to the 


eigenvalue decomposition (EVD) at less computational cost is shown. 
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IV. THE CRV DECOMPOSITION 


The purpose of this chapter is to describe the “Cross-product RV decomposition” 
[9]. The cross-product RV or CRV is a relatively new subspace tracking algorithm based 
on the URV decomposition [29]. In various signal processing applications, but mostly in 
noise suppression techniques, it is often needed to separate the signal subspace from the 
noise subspace. The most powerful tool to perform this task, as far as computational ac- 
curacy is concerned, is the singular value decomposition (SVD). Although the SVD is a 
very reliable means of determining the numerical rank of a correlation matrix, it has an 
important drawback. SVD is computationally expensive, especially for large matrices, 
and it is nonrecursive. This significant disadvantage can be effectively addressed by the 
CRV. 
A. CRV DECOMPOSITION 


We start by considering an Nx M (N=M) data matrix X of numerical rank r. 
Here, numerical rank means the number of singular values that are larger than a certain 


threshold. Therefore if o,,0,,......,0,, are the singular values of X, then 





O,20,2.....206, 2 threshold >> 0, % 0...) ®....% Oy ¥0. (4.1) 
By definition the singular values o, are directly related to the eigenvalues 4, of the data 


correlation matrix X¥"X since A, =o;. Thus the number of eigenvalues above the speci- 
fied threshold is the numerical rank of the matrix. In signal processing if the gap at the 
threshold is well determined, we say that the r largest eigenvalues are associated to the 
signal subspace and the M@—r smaller eigenvalues are associated to the noise subspace. 
In this way, each subspace is spanned by the corresponding eigenvectors. The CRV de- 
composition as already mentioned is based upon the URV algorithm [29], which is a 
rank-revealing decomposition proposed as an alternative to the SVD. The URV decom- 
position can provide a reliable estimate of the desired subspaces at a reduced computa- 
tional cost and may be written as 


ree ee oe 
X =URV" =U ye (4.2) 
0 G 
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where the matrices U, V are unitary, the matrices R, E, G are upper triangular, and F is 
a matrix of very small numbers. The smallest singular value (or eigenvalue) of the E ma- 


trix is oand the signal subspace is concentrated in this matrix. The noise subspace is 


concentrated mainly in the G matrix. 


The CRV decomposition is formulated by taking the product of the URV and its 
transpose. If as before X is an N x M data matrix, then [9] 


@ =X" xX =(URV") (URV")=VR" (U"U)RV". (4.3) 
Considering that U is a unitary matrix (U OP at ) and by setting 
W =R"R, (4.4) 


Equation (4.3) can be written as 


D=Vwv". (4.5) 
Equation (4.5) is a cross-product RV decomposition. By combining (4.2) and (4.4) the W 


matrix is written as 


W=R"R=|" FE F = ee ea = ee OF . (4.6) 
0 G||0 G F" G"||0 G F"E F"F+G"G 


Finally if d= E"E, B=E"F and C=F"F+G"G the rank-revealing CRV decompo- 





sition is formulated 
A B 
D=VWV" =VI ye. (4.7) 
BY“ 
Since the smallest singular value of E is o,, the smallest eigenvalue of A=E"E is 
A, =o. . Therefore the signal subspace power will be concentrated mainly in matrix 4 


and the noise subspace in matrix C. 
B. CRV UPDATING ALGORITHM 

In many signal processing applications, the rows of the data matrix X correspond 
to samples of a received signal corrupted with noise and the columns to different signals 


in time. If the signal and noise subspaces have already been estimated and a new signal is 
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received, then it would be more convenient to use the previously estimated noise power 
to approximate the new subspaces than to compute everything from the beginning. This 
procedure is called updating [29] and can be implemented recursively. The CRV updat- 


ing algorithm described next is a simplified version of the algorithm proposed in [9]. 
We start our description by assuming that the decomposition ® =VWV" is known 
and the updated one ® needs to be computed. The updated decomposition is given by 
D=VWV" =ab+ xx", (4.8) 
where @ is the fading (or forgetting) factor with values between 0 and | and x the newly 
received M x1 data vector. In our implementation, the fading factor is taken to be 1. Thus 
it can be omitted for the rest of the algorithm description. The next step is to define the 
intermediate vector z given by 
z=V" x, (4.9) 


so that equation (4.8) may be written as, 
@ =VWV" +Vz(Vz)" =VWV" +V22"V" =V(W+22")\V" =VEV", (4.10) 
and the problem is transformed into updating the matrix Y¥=W+2zz" . If the numerical 


rank of W according to a user specified threshold is r then ¥ is of the form (example for 


FS3) 




















x X% xe 2 e x 
XX xX €@ e x 
ot ee e x 
YW=l|e e eee e|+|x [x X X NX see ve x|; (4.11) 
eeeee e : ZH 
eeeee e x 
L alt ae 
W 2: 





where the upper left rxr elements of W, denoted with x’s, will have very large values 
compared to the rest (M—r)x(M—r) of the elements, denoted with e’s. Now that the 
matrix ¥ is defined, the next step in the algorithm is to find a Householder matrix P and 


apply a transformation to keep the first +1 elements of z by clearing the last M—r-1 
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elements. Householder matrices P (or elementary reflectors) [31] have the important abil- 


ity to zero specified entries in a matrix or vector and are defined as 


Pal-2-— (4.12) 
viv 
with J the identity matrix, 
v=zelel, an (4.13) 
and 
1 
0 
eZ=|. |, (4.14) 
0 
so that 
Pz=+'(lel, ae (4.15) 


For example if z= [2, 1, 4, 2] then zl, =5,v= [Fs 1, 4, 2) and 


[2 1 4 2 
L. -=34 7 - APT 2/7 
S44 AAT: SOFT. BT 7 
2. Died Sie Ssh 


> 





such that 


Pz=—5e, =[-5, 0, 0, OJ’. 


After computing the (M—r)x(M-—r) Householder matrix P, needed to clear the last 


M -r-1 entries of vector z, we form the M x M block diagonal matrix P as 


~ i 0,(M-r) 
P= (4.16) 
0, P. 


M-r)xr (M-r)x(M-r) 
With the help of equation (4.16) the updated matrix Y is given by 


Y=PY P= PW act | PX (4.17) 
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and is of the form 


at 
I 


(4.18) 


- RO SB BR & & 
/ DR SB BR BR 
“QR SB BR BR 
“QR BSB BR BR & 
vn 8 WR WR MN 
vn 8 WR WR MN 








ee ee a ae ee ata 
The upper left (r + 1)x (r +1) submatrix of ¥ denoted with x’s is unknown and, in order 


to proceed with the next cycle of the iterative algorithm, the rank 7 of ¥ needs to be es- 


timated. This is done by computing the eigenvalues of the (r+1)x(r+1) block of x’s 


and comparing the two smallest ones by putting a simple decision threshold. Specifically, 


if A,,, 20.14, then we determine that the rank has changed and 7 =r +1. On the other 


hand if 2,,,<0.1/2. the rank of the updated matrix Y remains unchanged and 7 =r. 


+1 
The way that the rank changes are determined in this implementation of the CRV updat- 
ing algorithm is not the most efficient one as far as computations are concerned. For “7” 
small (say 2 or 3) there are closed formulas that can be applied for the computation of the 
eigenvalues, as the Lanczos bidiagonilization algorithm [10]. These formulas offer the 
same accuracy at less computational cost but are not examined in this research. All the 
steps described above are repeated each time a new signal (vector x ) is received. At the 
end of each cycle, the correlation matrix ® is computed from Equation (4.10) and is 


given by 
D=VPV" (4.19) 
where 
V =VP. (4.20) 


A graphical representation of the CRV updating algorithm is illustrated in the fol- 
lowing flowchart (Figure 24). 
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Figure 24. CRV Updating Algorithm. 
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C. NUMERICAL EXAMPLE 

In this section a simple numerical example (done with MATLAB) is presented to 
demonstrate the accuracy of the correlation matrix estimate produced with the CRV up- 
dating algorithm illustrated in Figure 24. We start by creating a 10x10 data matrix X, 


Each column of X represents a random signal x,,i=1,2,....,10 of ones and minus ones 


corrupted with additive white noise: 








[-0.315 -0.781 -2.837 1.506 -2.731 -0.006 -0.223 0.654 -0.199 -1.379 
0.698 1.66 -0.502 -2.028 1478 0.306 -0.696 1.986 1.88 —1.242 
-1.599 2351 -0.18 --0.9 0.552 2.199 -1.73 1.643 0.775 —0.165 
0.219 1.136 0.765 0.833 1.386 -3.577 -2.143 1.919 -0.703 1.756 
ye 1.602 -0.833 -0.631 1.689 1.053 -1.086 -0.413 -2.248 -1.521 2.164 (421) 
0.057 -—0.705 -1.317 0.883 0.513 —-0.613 0.408 -0.842 0.844 —-2.023 
0.023 1.276 0.203 -0.674 1.244 0.138 1.518 1.788 0.901 0.701 
0.932 1.394 -0.309 0.904 1.718 2.23 2492 0.422 —-0.002 —-1.494 
1.081 -1.098 0.957 1.031 -0.846 1.641 0.913 -0.472 -0.565 —-0.827 
[—2.767 1176 1.324-° =1.613 --0.866> 0.095. 0.987. -2.671-- 0974. / 1354 


pF Xi My des aa sah ee die Ae bar 

















The correlation matrix estimate is computed for the signals x, and x, by going 


through all the steps of the first two cycles of the algorithm. 

if First Cycle 

The algorithm is initialized by setting W=0 andV =/. Both matrices are of di- 
mension 10x10 and the rank r of W is zero. Next we need to compute the z vector given 


by Equation (4.9): 


| 0.315 
0.698 
—1.599 
~0.219 


GEV GET eS ee : (4.22) 
—0.057 
—0.023 
0.932 
1.081 


| -2.767 
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where x, is the first column of the data matrix _X. From Equation (4.16) and the fact that 


r=0 we get P=P, where P is the Householder matrix needed to clear (set to zero) all 


but the first element of z. The Householder matrix P is computed with the help of the ex- 


pressions (4.12) through (4.14). Now that P is defined the updated matrix Y given by 


(4.17) is 


and the updated matrix V by (4.20): 


| 0.08 
0.177 
~0.406 
~0.055 
0.407 
~0.014 
~0.006 
0.237 
0.275 
| -0.703 


— 
Il 





0.177 
0.97 
0.066 
0.009 
—0.067 
0.002 
0.001 
—0.038 
—0.045 
0.115 


—0.406 
0.066 
0.846 

—0.021 
0.153 

—0.005 

—0.002 
0.089 
0.103 

0.264 


—0.055 
0.009 
—0.021 
0.997 
0.021 
—0.0007 
—0.0003 
0.012 
0.014 
—0.036 


0.407 
—0.067 
0.153 
0.021 
0.846 
0.005 
0.002 
—0.089 
—0.103 
0.265 


—0.014 
0.002 
—0.005 
—0.0007 
0.005 
0.999 
0.000 
0.003 
0.003 
—0.009 


—0.006 
0.001 
—0.002 
—0.0003 
0.002 
0.000 
0.999 
0.001 
0.001 
—0.003 


0.237 
—0.038 
0.089 
0.012 
—0.089 
0.003 
0.001 
0.947 
—0.06 
0.154 


0.275 
—0.045 
0.103 
0.014 
—0.103 
0.003 
0.001 
—0.06 
0.929 
0.179 


The last step is to compute the correlation matrix estimate @ (4.19): 


r 0.099 
0.22 
0.505 
0.069 
~0.506 
0.018 
0.007 
~0.294 
~0.341 
| 0.873 


Br 
lI 





0.22 
0.487 
1.117 
—0.153 
1.119 
—0.039 
—0.016 
0.651 
0.755 
1,932 


0.505 
—1.117 
2.558 
0.351 
—2.564 
0.091 
0.038 
—1.491 
1573 
4.426 


0.069 
—0.153 
0.351 
0.048 
—).352 
0.012 
0.005 
—0.205 
0.237 
0.608 


—0.506 
Td9 
2.564 
0.352 
2.569 
—0.091 
—0.038 
1.494 
1.734 
—4.435 
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0.018 
—0.039 
0.091 
0.012 
—0.091 
0.003 
0.001 
—0.053 
—0.061 
0.158 


0.007 
—0.016 
0.038 
0.005 
—0.038 
0.001 
0.0005 
—0.022 
—0.025 
0.066 


—0.295 
0.651 
—1.491 
—0.205 
1.494 
—0.053 
—0.022 
0.868 
1.008 
—2.579 


0.341 
0.755 
—1.73 

—0.237 
1.734 
—0.061 

—0.025 
1.008 

1.17 

2.993 


0 


10x10 


~0.703 | 
0.115 
~0.264 
~0.036 
0.265 
~0.009 
~0,003 
0.154 
0.179 
0.541 | 





0.873 | 
~1.932 
4.426 
0.608 
4,435 
0.158 
0.066 
-2.579 
~2.993 
7.656 | 





(4.23) 


(4.24) 


(4.25) 


The rank r of ¥ is determined (only in the first cycle of the algorithm) with the Matlab 


function “rank,” which provides an estimate of the number of linearly independent rows 


or columns of a matrix. In this example the rank r of ¥ is one. 


Z; Second Cycle 


The CRV updating algorithm is a recursive one; thus the updated matrices V and 


W =¥ computed in the first cycle are used to initialize the second cycle by setting 


W =W,V =V and r=1. The z vector is computed from (4.9): 


vou 
z=V"x,= 


| —2.388 | 
—1.396 
1.746 
1.053 
—0.226 
—0.726 
1.267 
1.747 
—0.689 








| 0.129 | 


(4.26) 


where x, is the second column of X. The Householder matrix P needed to zero all except 


the first two(r +1) elements of z is determined from Equations (4.12) through (4.14). For 


computational efficiency P is determined by operating on the last 9 (10 — r) elements of z 


and is of the dimension 9x9: 


[-0.405 0.507 0.306 
0.507 0.816 0.110 
0.306 -0.110 0.933 

0.065 0.023 0.014 

P=|-0.211 0.076 0.046 
0.368 —0.132 —0.080 
0.507 -0.183 —0.110 

0.200 0.072 0.043 

| 0.037 -0.013 —0.008 





—0.065 
0.023 
0.014 
0.996 

—0.009 
0.017 
0.023 

—).009 
0.001 


0.211 
0.076 
0.046 

—0.009 
0.968 
0.055 
0.076 

—0.030 
0.005 


oo 


0.368 
—0.132 
—0.080 

0.017 

0.055 

0.903 
0.133 

0.052 
—0.009 


0.507 —0.200 
0.183 0.072 
0.110 0.043 

0.023 —0.009 

0.076 —0.030 
0.133 0.052 

0.816 0.072 

0.072 0.971 
—0.013 0.005 


0.037 
~0.013 
0.008 
0.001 
0.005 
0.009 
0.013 
0.005 
0.998 |, , 


(4.27) 


The block diagonal matrix P is determined from Equation (4.16) as 


~ if On -r I 0 
P- (M ) -| 1 =| (4.28) 
0, F One Py 


M-r)xr (M-r)x(M-r) 


where P is the Householder matrix given by (4.27). The resulting matrix is 


0 0 0 0 0 0 0 0 0 

—0.405 0.507 0.306 —-0.065 -0.211 0.368 0.507 -0.200 0.037 
0.507 0.816 -0.110 0.023 0.076 -0.132 -0.183 0.072 —-0.013 
0.306 -0.110 0.933 0.014 0.046 —-0.080 -0.110 0.043 —0.008 
—0.065 0.023 0.014 0.996 —-0.009 0.017 0.023 -0.009 0.001 (4.29) 
0.211 0.076 0.046 -0.009 0.968 0.055 0.076 -0.030 0.005 |’ 
0.368 -0.132 -0.080 0.017 0.055 0.903 -0.133 0.052 —-0.009 
0.507 -0.183 -0.110 0.023 0.076 -0.133 0.816 0.072 —-0.013 
—0.200 0.072 0.043 —-0.009 -0.030 0.052 0.072 0.971 0.005 
0.037 -0.013 —0.008 0.001 0.005 -0.009 -0.013 0.005 0.998 | 


coo ccoocoo oF 








— 
oO 


and the updated matrix ¥ = PyP” = P(W +27" VeP produces the following: 


[ 21.17 -8.219 0 + 0 
—8.219 11.838 0 -: 0 
0 0 0 + 0 ; (4.30) 


“Er 
Il 








L *S Thoxio 
where the entries of ¥ , denoted with zeros, are very small numbers (< 10°) . The eigen- 


values of the upper left 2x2 (r+1xr+1) block of ¥ are 4,=25.955 and 4, =7.052 


and, sinceA,>0.1/,, we determine that the rank has changed by one and is 
nowr=r+1=2. It must be noted here that the decision threshold depends on the appli- 
cation, and there is no specific method known so far that could lead to an optimum choice 
[9]. The signal subspace power is concentrated in the upper left 2x2 block of Y whereas 


the noise subspace is in the lower right 8x8 block. Finally, the correlation matrix esti- 


mate ® calculated from Equations (4.19) and (4.20) is the following: 
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| 0.710 1.076 1.331 0.818 0.144 0.569 -0.989 -1.383 0.516 -0.045| 
1.076 3.244 -5.021 -—2.040 2.503 1.131 2.136 —1.664 2.579 —3.886 
—1.331 -5.021 8.087 3.023 -4.523 -1.566 3.039 1.788 -4.313 7.192 
0.818 2.040 3.023 1.339 —1.299 —0.788 1455 1.379 —-1.486 1.945 
0.144 2.503 -4.523 -1.299 3.263 0495 —-1.102 0.332 2.649 -5.415 (4.31) 
0.569 1.131 -1.566 -0.788 0.495 0.500 —-0.898 -1.036 0.712 —-0.671 
—0.989 -2.136 3.039 1455 —-1.102 -0.898 1.630 1.757 -1.428 1.567 
—1.383 —1.664 1.788 1.379 0.332 —-1.036 1.757 2.813 —0.523 —0.938 
0.516 2.579 -4.313 -1486 2.649 0.712 -1.428 -0.523 2.377 4.285 
| 0.045 -3.886 7.192 1.945 -5415 -0.671 1.567 -0.938 -4.285 9.040 | 


Ba 
I 








The exact correlation matrices for the first two signals x,and x,are determined 


with the following expressions: 


R= sa, (4.32) 


and 
RS R45, (4.33) 


After carrying out the calculations, the results from Equations (4.32) and (4.33) 
are exactly the same with the correlation matrix estimates (4.25) and (4.31) computed 
with the CRV updating algorithm. The eigenvalues of the exact correlation matrix and the 
eigenvalues of the estimated one for the 1“, 3™ 6" and 10" iteration are illustrated in 
Figure 25. By looking at the figure, we conclude that the estimate produced by the CRV 


algorithm is very accurate. 
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Eigenvalues of Correlation Matrix Estimate Computed With the CRV Algorithm 
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Figure 25. Values of the m-th Eigenvalue of the Estimated and the Exact Correlation 
Matrix. 


In this chapter, the CRV updating algorithm [9] for subspace tracking was de- 
scribed and a simple numerical example was presented so that its accuracy could be 
evaluated. In the next chapter of this thesis the performance of the algorithm is examined 
by implementing it in a previously proposed communication scheme [2] under varying 


channel conditions. 
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V. IMPLEMENTATION AND PERFORMANCE OF THE CRV 
UPDATING ALGORITHM 


This chapter describes the implementation of the CRV updating algorithm in a 
previously proposed direct sequence spread-spectrum (DS-SS) communication scheme 
[2]. The performance of the algorithm is evaluated under different channel conditions and 
design parameters (packet size, bit rate). The theory relevant to spread-spectrum (SS) 
signaling is presented first and the receiver structure is described next. The performance 


results are presented last in terms of Bit-Error-Rate (BER) vs. Bit Energy per Noise 


Power (£,/N,,), as it is usually done in all communication systems. 


A. DIRECT SEQUENCE SPREAD SPECTRUM (DS-SS) 

Spread-spectrum (SS) techniques are preferred in military communications be- 
cause of their low probability of intercept (LPI), low probability of detection (LPD) and 
their anti-jam capability [8]. In SS applications the resulting bandwidth of the transmitted 
signal is much greater than the information bandwidth and the transmitted signal appears 
to an unauthorized receiver as noise. The two most commonly used SS signaling schemes 
are frequency hopping (FH) and direct sequence (DS). Only DS is described in this the- 
sis. The block diagram of a simple DS-SS system is illustrated in Figure 26. 


ee 
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Figure 26. A Simple DS-SS System. 


The binary message signal is spread at baseband by a locally generated pseudo- 


noise code (PN-code). The PN sequence c(t), also known as the chipping sequence, is a 
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binary signal (c(t) = +1) that appears random and is produced at much higher frequency 


than the data that is to be transmitted. That is, 


R.>>R,, (4.34) 


where R, is the chip rate in chips-per-second [cps] and R, is the information data rate in 
bits-per-second [bps]. Since R, >> R, the chip duration 7. is much less than the bit dura- 
tion 7, . The processing gain (or multiplicity factor) N (always an integer) in chips-per-bit 


[cpb] for a DS system is 





Naa, (4.35) 


Therefore the message spectrum is spread N times before the message is modulated and 
transmitted. The wide bandwidth provided by the PN code allows the signal to be hidden 
in noise since its power is dropped below the noise threshold. In this way SS signals are 
very difficult to detect (LPD). The process of spreading by multiplication with a PN se- 
quence is shown in Figure 27 for the case of N= 2 (T, =0.57,). 
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Figure 27. Example of Spreading in a DS-SS System. 
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There are many different PN codes used in DS-SS communications like the 
maximal length codes, the Hadamard-Walsh codes, the Gold codes, the Kasami codes, 
and others. In the implementation of this thesis research, the spreading of the data signal 
is done by Gold codes, which are produced with the Gold code generators of 
SIMULINK. Gold codes have very good cross-correlation properties (the cross- 
correlation of two different Gold codes is very small) and are preferred when multiple 
users must share the same frequency spectrum (each user is assigned a different code). As 
shown in Figure 21, the exact replica of the PN code used at the transmitter is needed in 
order for the receiver to decode the transmitted signal. The decoding is done by a simple 
multiplication (de-spreading) of the received signal with the chipping code 


sincec’(t) =1 (c(t) = =a A The fact that the receiver must be able to generate the same 


spreading sequences in a deterministic way, just as the transmitter does, provides a cer- 
tain degree of communication security. The transmitted signal if detected by a non- 
intended receiver cannot be intercepted (decoded) without prior knowledge of the spe- 


cific PN sequence used for the spreading of the information signal (LPI). 


The third important advantage in using SS signaling is the ability to resist hostile 
interference (jamming). A useful parameter in specifying the performance of a spread- 
spectrum system in the presence of jamming is the processing gain N. Since the message 
signal is spread before transmission, the interference bandwidth occupies only a small 
portion of the SS signal bandwidth, as shown in Figure 28 (left plot). At the receiver, af- 
ter de-spreading the received signal by multiplying with the same PN code (the one used 
at the transmitter), the spectrum of Figure 28 (right plot) is produced. Most of the origi- 
nal interference energy is eliminated by spreading (except for the energy corresponding 
to the interference spectrum overlapping with the signal spectrum) and the original signal 
can be recovered. It is also clear from the same figure that the greater the processing gain 
N of the SS system, the better is the performance of the system against interference (bar- 


rage noise jamming) [18]. 
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Figure 28. Spectrum of the Transmitted SS Signal (left plot) and of the Received Sig- 

nal (right plot) after De-spreading [From Ref. 18]. 

B. PREVIOUSLY PROPOSED DS-SS SCHEME 

The channel equalization proposed by [2] is based on multi-rate signal processing 
and subspace decomposition. The analysis is done in the discrete time domain and the 
transmitted signal is assumed to be a binary sequence (+1, -1), which is then chipped by a 
short PN code and transmitted through the channel. Also, error correction coding and 
modulation were not examined and are not implemented in this thesis. Thus the 
worst-case scenario is simulated since both processes (error correction coding and modu- 
lation) can only improve the equalization mechanism. The block diagram of the simpli- 


fied DS-SS system is shown in Figure 29. 
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Figure 29. Simplified DS-SS System [After Ref. 2]. 


The information sequence m[n] is first up-sampled by N (add N—1 zeros be- 


tween bits) in order to match the chipping rate R, and then is spread by convolution with 
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the chipping sequence c[n]. Spreading can also be done by multiplying the up-sampled 
(where each bit is repeated N times) data sequence m[n] with c[n]. The spread signal 
s[n] is then convolved with the channel impulse response h[n] and noise w[njis added 


(white or colored) to account for the random interference. Subsequently, the received se- 


quence is de-spread at the receiver by convolution with a time-reversed replica c[—n] of 
the PN-code c[n] and then is down-sampled by WN. In order to effectively remove the 


““pseudo-randomization” imposed on the data sequence by the spreading code, a convolu- 


tion between c[n] and c[—n] must yield a delta function, 
c[n]* c[—n] = 6[n]. (4.36) 
If c[n] =[cg, ¢,, +5 Cy], where cy, ¢,, ....., Cy_; are the chips of the PN-code, then 
c[-n]= leis a or |. (4.37) 


For example, when a Gold code c[n] of length 125 is used and c[—n] is given by Equa- 
tion (4.37), the result of (4.36) is the impulse shown in Figure 30. 














Figure 30. The Result of c[m]*c[—n] When Gold Codes are Used. 
The equivalent state-space structure of the communication system shown in Fig- 


ure 29 is illustrated in Figure 31. 
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Figure 31. State-Space Structure of the DS-SS System [After Ref. 2]. 


The spreading of the information sequence m[n] is accomplished by filtering the 


up-sampled (by NV) data sequence by a N-th order linear filter defined (in the z-domain) as 
C(z) =e) +6, 2° 46,27 tut Oy, 20. (4.38) 


The spread signal s[n]is then filtered by a second N-th order linear filter (channel im- 


pulse response), 


H(z)=hythz'+hz?+....+hy z”, (4.39) 
to account for the corruption caused by the underwater acoustic channel. At the receiver 
the received sequence is de-spread by multiplication (and not convolution) with a replica 
of the PN-code used at the transmitter. In practice the length of the received sequence 
does not change after de-spreading. The summation step implied by the convolution ap- 


proach (c[n] * c{-n]) is performed later at the equalizer when the estimate m[n] of m[n] 


is computed. At the equalizer (see Figure 31), the de-spread sequence is forward delayed 
(z elements in the state-space structure) and down-sampled by N . The result of this proc- 
ess is N different versions of the original information sequence (corrupted by the fading 


channel) denoted asy,[n], y,[7],....... ¥y in]. For a better understanding of how the 


structure of Figure 31 works, simple examples are presented in Figure 32 (transmitter 


section) and Figure 33 (channel and receiver section). 
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Transmitter Section * 
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Figure 32. Up-sampling and Spreading by Multiplication with the PN-code. 


The data sequence m[n]consists of 5 bits [1, 1, —1, 1, —1] and the processing gain 
is assumed to be three( NV =3). Since N =3, mln] is up-sampled by 3 (each bit is re- 
peated three times) and then bit-by-bit multiplied with the PN-code c[n] to yield the 
spread sequence s[n]. Subsequently, white noise w[n] is added to account for the chan- 
nel corruption. The sequence s[”]+w[n] is then de-spread at the receiver by multiplica- 
tion with a replica of c[n] (Figure 33, third plot from top). Finally, the de-spread se- 
quence is forward delayed and down-sampled by WN, so that three different ver- 
sions y,[n], y,[”], y,[”] of the original data sequence are produced at the equalizer input. 
The three sequences y,[n], y,[”], y,[”] are shown in Figure 33 (first plot from bottom) in 


different colors (blue, red, and green respectively). 
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Channel and Receiver Section 
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Figure 33. Channel Corruption, De-spreading and Different Versions of the Original 
Data Sequence. 


The vectors y)[7], ¥[7],..... Vy["] (see Figure 31) form the data matrix 


yin] 


y [7] 


vnl=| . oil (4.40) 


Vyaln] 


which is then used to determine the signal and noise subspace (the most important part in 


the equalization process). The estimated sequence m[n] is given by 


m[n]= f" yn], (4.41) 
where f is the optimum filter that maximizes the signal-to-noise ratio at the receiver. 


The channel equalization is performed with two different blind equalization algorithms, 


depending on whether the synchronization between the receiver and the transmitter was 
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achieved or not. The two algorithms are based on matched filter theory and are briefly 
described in the following subsections. A detailed derivation and description can be 
found in Reference [2]. 

1. Demodulation with Synchronization 

This is the case where the replica of the PN-code used at the transmitter is com- 
pletely aligned (in time) with the received sequence. Thus the de-spreading is performed 


correctly. 


a. Multipath with Additive White Gaussian Noise (AWGN) 
The received sequence y[n] is equal to s[n]*A[n]+w[n] where s[n] is 
the spread by the PN-code information sequence, w[n] is the AWGN, and /[n] is the 


channel impulse response. It must be noted here that the maximum length of the impulse 


response /[n], according to the theoretical approach described in the previous section, is 
N, where N is the processing gain (see Figure 31). If m[n] is an estimate of the informa- 


tion sequence m[n] then 


m[n]= f" - y[n], (4.42) 
where f is the eigenvector corresponding to the maximum eigenvalue of the data corre- 


lation matrix R,,, given by 
R,, =E{y[nly" [n}}. (4.43) 


b. Multipath with Additive Colored Noise (ACN) 


The received sequence y[n] is equal to s[n]*h[n]+v[n] where v[n] is the 
additive colored noise. In order to be able to follow the approach mentioned above in 


(4.42), a simultaneous diagonalization of the correlation matrices of the signal R,, and 


the noise R,, is needed. This is done with help of the Mahalanobis transformation, de- 


fined as [32] 


y[nj=M yn], (4.44) 
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where M~"” is the transformation matrix given by 

MW? = BANE =R,”, (4.45) 
and E,, A, are the matrices produced from the eigenvalue decomposition of 
R,=E {vial [n]} . Specifically £, is a matrix with columns of the eigenvectors of R,, 
and A, is a diagonal matrix whose elements are the corresponding eigenvalues. The cor- 


relation matrix R\, of y'[n] (4.44) is computed from 


R',= E{ yin (y [ny } =E (wr y{n])(M"? vial)’ =M"?R M", (4.46) 
and, since M'* = Ro , from 


1 _ p-l/2 -1/2 
Ry, =R, RR, (4.47) 

Because of this whitening transformation in the new coordinate system the 
colored noise is considered white and the estimate m[n]of the data signal m[n] is given 


by (4.42) 


in| =(f") yn], (4.48) 
where f’ is the eigenvector corresponding to the maximum eigenvalue of the autocorre- 
lation matrix R’,. 

2 Demodulation without Synchronization 

The underwater communication channel is very challenging and synchronization 
is not always possible. The slightest misalignment in the PN-code can incorrectly de- 
spread the received sequence. When synchronization is achieved, the data correlation ma- 


trix R,, (4.43) has only one dominant eigenvalue and the signal subspace is easily dis- 
criminated from the noise subspace. However, this is not the case when synchronization 
is not possible. Suppose that, at a given time, y[”] contains two symbols, the current 


symbol m[n] and a portion of the previous symbol m[n —1], i.e. 


See RE eC al (4.49) 
— = = E g | m[n-1] 
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where g,, g, are vectors depending on the autocorrelation of the PN-code and the com- 


munication channel. In this case y[7] spans a two-dimensional subspace defined by the 


vectors g, and g, and there are two dominant eigenvalues instead of one in the data cor- 


relation matrix R,,. If A,[”], 2,[n] are the two dominant eigenvalues and 4, @ their re- 


spective eigenvectors, y[7] is given by: 


yn} = @Aln]+ 24 [7]. 
Since the eigenvectors are orthonormal (¢,"e, =1 and e,'e, =0), it is true that 


A[n] =e": y[n] 


and 


A,[n]= 2" y(n]. 


By combining Equations (4.51) and (4.52) in matrix form, we get 


Aln]}_je' | rel 
al = El yln] = A[n]= El yin]. 


Substituting (4.49) into (4.53) yields 


eer | mn] 
an" ][6 8}, 


2 


, mn] 
If we solve Equation (4.54) for 
mi[n— 


min] = e’ | ‘i min] 7 kel! 
Eel Ls s]| ae gl lar ae 


where we assumed that g,, g,(and &,’, k, ) are linearly independent, and 


El -(E} sl). 


From (4.55) the two information symbols are given by 


mn] =k; -Afn], 
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(4.50) 


(4.51) 


(4.52) 


(4.53) 


(4.54) 


(4.55) 


(4.56) 


(4.57) 


and 


m[n-l=k, -A[n] S nf{nj=k - An +1]. (4.58) 
A combination of (4.57) and (4.58) yields 


Ki -A[nl=ky -A[n+1] Ok! -Aln]-k -A[n+1] = 0, (4.59) 
and in matrix form 


» arf all]. 
ES Ea Pane (4.60) 


k, 
In order to calculate k -| = 


| we start by calculating the autocorrelation ma- 
mee) 


trix R,, = > 


n 


Aln] 


[2"[n] 4"[n+1]| with help of Equation (4.53). From (4.60) it is 
A{n+1] 


obvious that 

|e? -ke |-R,, =. (4.61) 
Therefore, k is equal to the eigenvector corresponding to the minimum eigenvalue of 
R,,. Once the estimate for k is known, the estimated sequence m[n] is derived from 
(4.55). The equalization algorithm mentioned above is tailored for the specific case that, 
at any given time, [7] contains only two symbols (slight misalignment). If this is not the 
case and y[n] contains more than two symbols, the same approach can be followed (in 


this case there are more than two dominant eigenvalues) that would lead to a similar algo- 
rithm. 
C. PROPOSED DS-SS SCHEME 

The combination of the DS-SS modulation scheme and the receiver structure 
based on the aforementioned blind equalization algorithms was proven very robust as part 
of a previous thesis [2]. However, there is a major problem associated with this scheme. 
The desired subspaces can only be estimated offline since the whole received packet is 


needed in order to calculate the data correlation matrix R,,. This places a computational 


burden that increases as the packet length increases. 
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In this thesis the signal and noise subspaces are estimated online (recursively) 
with help of the CRV updating algorithm presented in Chapter IV. By following this ap- 
proach, the desired subspaces are efficiently estimated in a recursive manner. This is 


done by feeding each column of the received data matrix y[n] (4.40) into the CRV algo- 
rithm. The columns of y[n] represent the received spread bits, corresponding to one 


information bit spread by the PN-code and corrupted by the fading channel. 


In the previously proposed communication scheme, the channel impulse response 
was generated using the Bel/hop underwater acoustic propagation model [5]. The time 
variability of the channel was introduced to the static impulse response by allowing each 
of the multipath components to vary randomly within an arbitrary range. The random 
phase shift caused by the micro-paths was approximated by adding a random phase factor 


to each of the amplitudes. 


In this thesis the channel impulse response is modeled using the Monterey-Miami 
parabolic equation model (MMPE) [7] described in Chapter III. The broadband impulse 
response is generated for the frequency range 11,500+1024Hz (BW = 2048Hz) under 


various perturbations including the influence of Doppler due to source motion. The 
worst-case scenario is examined by evaluating the CRV algorithm using a time-variable 
impulse response that also accounts for the interface roughness, the turbulence, the inter- 
nal waves, and the volume scattering commonly found in shallow water environments. 
D. SIMULATION RESULTS 

The performance of the CRV updating algorithm evaluated under different chan- 
nel conditions and design parameters (packet size, bit rate, etc.) is examined next. The 
DS-SS system structure of Figure 31 was implemented in MATLAB and the simulation 


results are presented in terms of Bit-Error-Rate (BER) vs. Bit Energy per Noise Power 
(E,/ N,). Each simulation involved sending an adequate number of packets for each 
value of E,/N, in order to create sufficient statistics. The E,/N, is a scaled version of 


the signal-to-noise ratio (SNR) defined as 


ot eRe (4.62) 


0 b 


69 


where BW is the available communication bandwidth and R, is the bit rate. The range of 
the desired £,/N, used for the various simulations was between 0 and 30 dB. The BER 
(probability of a bit error) is a metric used to describe the reliability of a digital commu- 


nication system quantitatively and is defined as 


ER= Total Number of sis . (4.63) 
Total Number of Bits Sent 


1. Performance in AWGN 


The received sequence is s[m|]+w[n] where s[n] is the spread information se- 
quence and w{n] is the AWGN. In this case the channel corruption is not taken into ac- 


count; thus, the results of this simulation can be treated as a benchmark for evaluating the 
receiver performance. Figure 34 shows the resulting BER for different packet sizes and 


Figure 35 the BER for various data rates R, . 
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Figure 34. . BER in an AWGN Channel When Varying the Packet Size 
(R, =40bps, R, = 2400cps). 
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Figure 35. BER inan AWGN Channel When Varying the Data Rate R, 
(packet size = 72 bits, R, =2400cps). 

From the results presented above two observations can be made. First, we see 
that, if longer transmission packets are used, the performance is not degraded in this ideal 
case. The same conclusion stands for the case where a higher data rate is used. The per- 
formance seems to be better for a packet of 144 bits (Figure 34) and for a data rate of 60 
bps (Figure 35). This is possibly because of statistical variations. 

2. Performance in Additive Color Noise 

The colored noise is produced by filtering white noise by the 7"-order FIR filter 
given by 


1 


Q(z) = z° —0.92° + 0.4z* —0.2z7 +0.1z? -0.1z 





(4.64) 


This is the same filter used in [2] so that the performance results may be compared. 


Figure 36 illustrates the performance of the new subspace decomposition algo- 
rithm (CRY) in additive colored noise when varying the packet size. The results of this 
simulation are worse than the white noise case by approximately 4 dB. Figure 37 demon- 


strates the receiver performance for a 72-bit packet when the bit rate is varied in the range 
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(R, =40bps, R, = 2400cps). 
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of 40 to 100 bps. It is shown once again that the algorithm is stable even for a bit rate of 


100 bps. 
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Figure 36. 
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Performance in Additive Color Noise When Varying the Data Rate 
(packet size = 72 bits, R, 


Figure 37. 


J Performance in Multipath and AWGN 

The received sequence is s[n]*A[n]+w{n] where s[n] is the spread by the 
PN-code data sequence, w[n] is the AWGN, and /[n] is the channel impulse response. 
The impulse responses used for the simulations are the ones computed with help of the 
MMPE model. The performance of the CRV algorithm was evaluated for all the different 


perturbations described in Chapter III and the results are shown in Figure 38 through Fig- 
ure 44. 


As already mentioned, the communication scheme of Figure 31 is tailored for a 
channel impulse response equal to the processing gain WN. For all the following simula- 


tions, R, =40bps and R, =2400cps, thus N is 60. This is the reason all results are pre- 


sented for a channel length of at least 60 time samples. 


Among the individual perturbations, the most challenging are those corresponding 
to the interface roughness (Figure 39) and turbulence (Figure 42). In those cases the re- 


sults are good only for a channel length of 60 time samples. 


When all the examined perturbations are combined (Figure 43) the results are 
poor and on the average worse than the previous cases by 18 dB (for an impulse response 
of 60 samples). The results of the worst-case scenario where the Doppler is also taken 
into account are presented in Figure 44. In this case, equalization cannot be performed, 
even for a channel length of 60 time samples. However, it must be noted that the exam- 
ined perturbations are expected to be realistic but strong perturbations. It is doubtful that 


all of them combined will be encountered in many shallow water environments. 
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Figure 38. 
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Interface Roughness Perturbation: BER in a Fading Channel with AWGN 


Figure 39. 
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Figure 40. 
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Figure 41. 
When Varying the Channel Length (packet size = 72 bits, 


2400cps). 
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Figure 42. Turbulence Perturbation: BER in a Fading Channel with AWGN When 
Varying the Channel Length (packet size = 72 bits, R, =40 bps, R, =2400 cps) ‘ 
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Figure 43. All Perturbations Combined (no Doppler): BER in a Fading Channel with 
AWGN When Varying the Channel Length 
(packet size = 72 bits, R, =40 bps, R, =2400 cps). 
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Figure 44. All Perturbations Combined (with Doppler): BER in a Fading Channel 
with AWGN When Varying the Channel Length 


(packet size = 72 bits, R, =40 bps, R, =2400 cps). 


4. Performance in Multipath and Additive Color Noise 

This is the same case as the previous one but colored noise was added instead of 
white. The equalization algorithm used in this case is described by Equations (5.11) 
through (5.15) where once again the CRV is used to estimate the signal and noise sub- 
spaces. The resulting BER for a channel length of 60 time samples is shown in Figure 45. 
The algorithm works well for all the individual perturbations but not when all of them are 


combined. 
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Figure 45. Performance in a Fading Channel with Additive Color Noise When Vary- 
ing the Channel Length (packet size = 72 bits, R, =40 bps, R, =2400cps). 


3: Performance When Synchronization is Lost 

When the PN-code in the transmitter and its replica in the receiver are synchro- 
nized (aligned in time), the received sequence is despread correctly. However due to the 
volatile underwater environment, this may not always be the case. When synchronization 
is lost, the dominant eigenvalues of the received covariance matrix are more than one. 


When at any given time y[n] contains two symbols (the current symbol m[n] and a por- 


tion of the previous symbol m[n—1]), Equations (4.49) through (4.61) may be used for 
the recovery of the transmitted sequence. Therefore, in order to simulate the loss of syn- 
chronization, a random delay (of duration up to half a bit) was added to the received sig- 
nal. The performance results of the algorithm are poor for all individual perturbations. 
For the case in which all perturbations are combined and subsequently the Doppler Effect 


is added the resulting BER’s are shown in Figure 46. 
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Figure 46. BER for the Case in which Synchronization of the PN-code is Lost. 


6. Summary of Performance 
In this chapter, the performance of the three different equalization algorithms, 
with the CRV used for the estimation of the signal and noise subspace, was examined un- 


der various design parameters and channel conditions. 


When synchronization between the spreading sequence and its replica is achieved, 
the results of the proposed communication scheme are good for all the individual 
perturbations. For a fading channel with AWGN, a BER of 10° is possible for 


E,/N, 2 13dB whereas in a fading channel with additive color noise the same BER is 


achieved for E,/N, = 17dB. 


Among the examined perturbations, the ones corresponding to the interface 
roughness and the turbulence have more impact on the performance of the communica- 
tion scheme than the others. The small-scale fluctuations of sound speed within the water 
column caused by turbulence, and the scattering of the acoustic energy interacting with 
the rough bottom create a very challenging environment for underwater wireless commu- 


nications due to the greater breakdown in the coherence of the propagating field. 
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The performance of the algorithms is not the desired one when all the perturba- 
tions are combined and Doppler is added to account for the source motion. But as already 
mentioned the modeled perturbations are rather strong, creating a shallow water envi- 


ronment more challenging than usual. 


When synchronization is not achieved and the algorithm described earlier is used 
for the recovery of the information sequence, the results are very poor. As shown in 


Figure 46, a BER < 10° cannot be achieved even for E,/N, = 25dB. In the last chapter, 


the important findings of this thesis research are summarized and possible areas for fol- 


low-on work are presented. 
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VI. CONCLUSION 


The primary purpose of this thesis was to evaluate the performance of the CRV 
updating algorithm developed in a previously proposed communication scheme intended 
for use in an underwater acoustic network operating in a shallow-water environment. The 
transmitter and receiver were implemented in MATLAB and the underwater channel was 
modeled using the Monterey-Miami Parabolic Equation model (MMPE). Various simula- 
tions were performed under different channel conditions (interface roughness, turbulence, 
etc.) and design parameters (packet size, bit rate). 

A. FINDINGS 

The CRV was proposed [9] as an alternative to the eigenvalue decomposition 
(EVD). The main advantage offered by the CRV is that it can be updated online at less 
computational cost compared to the EVD. As shown in Chapter IV, the estimates of the 
signal and noise subspace produced by the recursively updated algorithm are very accu- 


rate. 


In order to examine the performance of the CRV, three distinct cases were simu- 
lated; a shallow-water fading channel with AWGN, a fading channel with .colored noise, 
and the case in which synchronization between the receiver and the transmitter were not 


achieved. Under all individual perturbations described in Chapter II, the equalization 


algorithm performed well and achieved a BER of <10° for values of £,/N, = 13 dB 


and E,/N,=17dB for white and colored noise, respectively. 


Among the modeled perturbations, the ones corresponding to the interface rough- 
ness and turbulence were the most challenging for the proposed equalization algorithm. 
As mentioned in Chapter HI, these perturbations cause the greater breakdown in the co- 


herence of the propagating field resulting in a severe multipath environment. 


On the other hand, when all perturbations were combined and the Doppler effect 
due to source motion was considered, the resulting bit-error-rates were not acceptable. 
For the case of no synchronization between the spreading and de-spreading sequence, the 


performance of the acquisition algorithm was also very poor. 
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B. FUTURE WORK 

The channel impulse response was obtained by post processing the output of the 
MMPE acoustic propagation model, which represented the response of the simulated un- 
dersea environment. Therefore, the proposed communication scheme was tested under 
simulated conditions. It would be very interesting to evaluate its performance using real 
impulse responses obtained from experimental data in a variety of shallow-water chan- 


nels. 


In the current implementation, the desired subspaces (signal, noise) were esti- 
mated online with help of the recursively updated CRV algorithm. However, the equali- 
zation based on matched filter theory is still performed offline. More research on design- 
ing an effective communication scheme that could operate online is needed. This can be 
partially achieved by using training bits in each transmission or by using the first packet 


as a training transmission. 


The CRV updating algorithm could be even more computationally efficient if the 


Lanczos bidiagonalization algorithm were used to compute the eigenvalues of the 'P ma- 


trix given by Equation (4.18). 


There is a family of rank-revealing orthogonal decompositions, referred to as 
UTV decompositions [33], offered as alternatives to the SVD. These algorithms are able 
to provide all the necessary subspace information, and most importantly, can be updated 
online. A comparison of their performance, with the results obtained with the CRV de- 
composition, could lead to the design of a better signal processing algorithm for the re- 


covery of the information signal in a multipath channel. 


Finally, the performance of the proposed transmitter/receiver structure should be 
evaluated as part of a complete communication system with error correction coding and 


modulation. 
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APPENDIX. MATLAB CODE 


The MATLAB source code used for the simulation of the communication system 
structure of Figure 31 and the implementation of the CRV updating algorithm is pre- 
sented in this appendix. The overall controlling m-file is called CRV_simulations.m. All 
the other necessary functions (m-files) are called by running this file and are listed bel- 


low. 

CRV.m — This function implements the CRV updating algorithm described in 
Chapter IV. The output of this function is the recursively updated correlation matrix @ , 
and the matrices W and V needed for the initialization of the next step of the algorithm. 


up.m — This function is used to upsample the message signal to the chipping rate 


so that the PN-code may be appended. 


spredesp.m — This function is used to spread the binary message signal with the 
Gold code, or de-spread the received sequence with the replica of the Gold code used at 


the transmitter. 


power_of noise.m — This function calculates the noise power associated with a 
given E,/N, [dB], a bit length 7, [s], and a sampling period 7, [Hz]. The output is used 


to generate the white or colored noise corresponding to the noise power. 


house.m — This function computes the Householder matrix P required to clear 


specific values of the received sequence (a required step in the CRV updating algorithm). 
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OBES ISIC I IIIA ICICI CIR CIE 


% Overall Controlling Program 

% T/R SIMULATION 

% Note: all secondary variables and needed 

% functions are called by executing this m-file 
% 

% developed by Pavlos Angelopoulos Mar 2004 
% last modified 4/13 

O62 He 2s ae 2 6 2 2 fe 2 he 2 2 he 2 6 2 6 22 fe 2 he 2 2c he 2s fe 2s fe 2 2s aie 2 ie 
clear; 

cle; 

% Offer options for different simulations 
disp(‘Choose one of the following options:'); 
disp(‘1. White Noise only’); 

disp('2. Color Noise only’); 

disp('3. Fading ChanneltA WGN’); 

disp('4. Fading Channel+ACN’); 

disp('S. Loss of Receiver Synchronization’); 





option=input(' '); 

disp("'); 

if option==3|option==4; 
% Offer options for different perturbations 
disp(‘Choose one of the following perturbations:'); 
disp(‘1. Basic’); 
disp('2. Interface Roughness’); 
disp(‘3. Internal Wave’); 
disp('4. Volume'); 
disp('5. Turbulence’); 
disp('6. All Perturbations (without Doppler)'); 
disp('7. All Perturbations (with Doppler)'); 
pert=input(' '); 
disp(‘'); 
% Load the desired Impulse Response 
if  pert==1; load ImpRes_basic; 
elseif pert==2; load ImpRes_intrough; 
elseif pert==3; load ImpRes_intwave; 
elseif pert==4; load ImpRes_volume; 
elseif pert==5; load ImpRes_turbulence; 
elseif pert==6; load ImpRes_allpert; 
elseif pert==7; load ImpRes_ doppler; 
end % end of "if pert" 
disp(*'); 

end % end of "if option==3|option==4" 





if option==3|option==4|option==5; 
disp('Choose the Channel Length:’); 
channel_length=input(' '); 
disp(‘'); 

end % end of "if option==3|option==4|option==5" 











disp('‘Choose # of Simulations:’); 
trials=input(' '); 
disp("'); 
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Rb = 40; % bit rate [bps] 

Re = 2400; % chip rate [cps] 

BW=2*Re; % bandwidth of spread signal 

Tb = 1/Rb; % bit Period [s] 

Tc = I/Re; % chip Period [s] 

N = Re/Rb; % processing gain 

P=Tb/Tc; % Upsampling required to match chip sequence 

load pn_code; % load PN chipping sequence "pn_code.mat" 
packet_length = 72; % number of bits transmitted per packet 
number_of_packets=1; % number of packets transmitted 


EbNo_dB=linspace(1,25,12); % determine Eb/No values to run simulations 
disp(‘START'); 


disp("'); 

Oot LO aS I POS a NA ee te A ee ee ee 
em oe IE ee ee eae eT 
% TRANSMITTER 

OS oc I re ate ao 2 as aR ed et ee as ene 
0) Re eS oar en re ee eee a SE 
for simulation=1 :trials; % #of simulations 

txt=sprintf(‘Simulation #%3d',simulation); 

disp(txt); 

tic; % starts clock to measure simulation run time 
txt2=sprintf(‘Simulation #%3d',simulation); 

h=waitbar(0,txt2); 


for n=1:length(EbNo_dB); = % run simulation for each of the EbNo values 
waitbar(n/length(EbNo_dB),h); 


upsampled_data_seq=up(data_seq,P); 
m=length(upsampled_data_seq); 


spreaded_signal=conv(upsampled_data_seq,pn_code(1:P)); 
spreaded_signal=spreaded_signal(1:m); 


% Calculate the noise power associated 
% with the simulated EbNo values 


noise _power=power_of noise(spreaded_signal,Tb,1/BW,EbNo_dB(n)); 
noise_std = sqrt(noise_power); % noise standard deviation 
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if option==1; 

white_noise=randn(1,length(spreaded_signal)); 

received _seq= spreaded_signal + white _noise*noise_std;% Add White Noise 
elseif option==2; 

white_noise=randn(1,m); 

color_noise=filter(1,[1 -0.9 0.4 -0.2 0.1 -0.1],white_noise);% filter white noise to produce color noise 
received_seq= spreaded_signal + color_noise*noise_std;% Add Color Noise 
elseif option==3; 
spreaded_signal=conv(imp_res(1:channel_length),spreaded_ signal); 
white_noise=randn(1,length(spreaded_signal)); 

received _seq= spreaded_signal + white _noise*noise_std;% Add White Noise 
elseif option==4; 
spreaded_signal=conv(imp_res(1:channel_length),spreaded_ signal); 
white_noise=randn(1,length(spreaded_signal)); 

color_noise=filter(1,[1 -0.9 0.4 -0.2 0.1 -0.1],white_noise); 

received_seq= spreaded_signal + color_noise*noise_std;% Add White Noise 
else 

load ImpRes_ doppler; 
spreaded_signal=conv(imp_res(1:channel_length),spreaded_signal); 
white_noise=randn(1,length(spreaded_signal)); 

received _seq= spreaded_signal + white _noise*noise_std;% Add White Noise 
end % end of "if option==1" 


if option==5; 

% Add random time t0 uniformly distributed between 0 and P 
% to simulate the Loss of Synchronization between the codes 
t0=round(rand*P/4)+1; 
received_seq=received_seq(t0:length(received_seq)); 

end 


despreaded_signal=spredesp(received_seq(1:m),pn_code(1:P)); 
despreaded_signal=despreaded_signal(1:m); 


Y=reshape(despreaded_signal,P,packet_length); 
columns=size(Y,2); 
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% Initialization 
W=zeros(P,P); 
V=eye(P); 
r=0; % Rank of matrix W (1=0) 
for k=1:columns; 
[phi_crv,W,V,rJ=CRV(Y(:,k),W,V,1); % call the "CRV.m" function 
end 


if option=1|option==3; 
[EigVec,EigVal] = eig(phi_crv); 
Opt_filter=EigVec(:,1); % Eigenvector of the largest eigenvalue 

elseif option==2|option==4; 
corrmatrix_color_noise=color_noise*color_noise'; 
[EigVec_color,EigVal_color]=eig(corrmatrix_color_noise); 
R=inv(sqrtm(EigVal_color)); 
R_sq_inv=EigVec_color*R*EigVec_color'; 
R_mah=R_sq_inv*phi_crv*R_sq_inv; 
[EigVec_mah,EigVal_mah]=eig(R_mah); 
Opt_filter=R_sq_inv*EigVec_mah(:,1); 

elseif option==5; 
[EigVec,EigVal] = eig(phi_crv); 





Opt_filter=EigVec(:,1); % Eigenvector of the largest eigenvalue 
EigVal=diag(EigVal); 
if EigVal(2)>0.2*EigVal(1); % Assume there is no synchronization 


lamda=EigVec(:,1:2)'*Y; 
E=[lamda(:, | :length(lamda)-1); lamda(:,2:length(lamda))]; 
R_lamda=E*E'; 
[Ve, De]=eig(R_lamda); 
lam_e=diag(De); 
I=find(lam_e==min(abs(lam_e))); 
ke=Ve(:,D; 
k1=ke(1:2); 
k2=-ke(3:4); 
end 

end 


if option==5 & EigVal(2)>0.1*EigVal(1); 
decoded_seq=sign([k1',k2']*E); 
else 
decoded_seq=sign(Opt_filter'*Y) ; 
end; 
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if option==5 & EigVal(2)>0.1*EigVal(1); 
diff=length(decoded_seq)-length(data_seq); 
if diff>0 
err=length(find(data_seq-decoded_seq(diff:length(data_seq)+diff-1)~=0)); 
elseif diff<0 
diff=abs(diff); 
err=length(find(data_seq(diff:length(decoded_seq)+diff-1)-decoded_seq~=0)); 
else 
err=sum(abs(data_seq-decoded_seq))/2; % # of errors 
end % end of "if diff>0" 
else 
err=sum(abs(data_seq-decoded_seq))/2;  % # of errors 
end % end of "if option==5 & EigVal(2)>0.1*EigVal(1)" 
if err>packet_length/2; 
channel_errors(n)=packet_length-err; 
else 
channel_errors(n)=err; 
end % end of "if err>packet_length/2" 


end % end of n loop 

close(h); 

toc; 
summary(simulation,:)=channel_errors./packet_length; 


end % end of simulations loop 
avg _errors=mean(summary); 


result=[EbNo_ dB; avg _ errors]; 


result(2,find(result(2,:)==0))=10%-6; 

figure; 
semilogy(result(1,:),result(2,:),'-',"MarkerSize',8); 

grid on; 

ylim([10%-5 1]); xlim([1 30]); 

ylabel('BER', 'fontsize', 13); 

xlabel('Eb/No (dB)','fontsize’, 13); 
legend(sprintf(‘Packet size: %3.0f bits',packet_length)); 


disp(‘STOP'); 

% offer option to save the simulation results 

% for post processing (plots) 

savedata=input('Save Simulation Results? (y or n)','s'); 

disp("'); 

if savedata == 'y'| savedata == 'Y' 

fileout=input('Enter Output filename (without extension): ','s'); 
eval(['save ' fileout ' result trials packet_length Rce;']); 


end 
if savedata == 'n' | savedata == 'N' 
end 
break 
o------------------------- END of Code------------------------------ 


function [phi_crv,W,V,r]=CR V(signal, W,V,r) 


OEE EO SSSI CGISCI CISC CICCICICR ICI CII ICICI 


% CRV UPDATED ALGORITHM 

% "signal" is the new signal used to update the correlation matrix 
% "W" and "V" are the matrices needed to initialize the algorithm 
% "r" is the rank of W 

% "phi_crv" is the esimate of the correlation matrix 

% 

© developed by Pavlos Angelopoulos Feb 2004 

0 last modified 3/28 


OEE SORES SCC CCI GICC ICICI ICICI ICI ICICI 


S38 


a=1; % Fading Factor 
P=length(signal); 
W=a*W; 

z=V"* signal; 


P_hat=house(z(r+1:end)); % Householder vector: keep r+1 element and zero the rest 
P_diag=[eye(r) zeros(r,P-r);zeros(P-r,r) P_hat]; % P=diag(Ir,P_hat) 
W=P_diag'*(W+z*z')*P_ diag; % update matrix W 
V=V*P_diag; % update matrix V 
warning off; 
eigenv=abs(eigs(W(1:r+1,1:r+1),r+1)); % compute eigenvalues of (r+1)x(r+1) block of W 
if eigenv(r+1)>0.1*eigenv(r); % Threshold to check rank change 
r=rt+1; 
end 
end 
phi_crv=V*W*V'; 
else 
W=W3z*z'; 
phi_crv=V*W*V'; 
r=rank(W); 
end 


function y=up(x,M) 


ORES O SIGS ICSI ISCIC CIC ICICI CIC 


% Upsample: A function to upsample a signal x by an integer M 
% 

% developed by Pavlos Angelopoulos Nov 2003 
% last modified 3/28 


ORES O RSIS CICCIGICICCICCIICICGIICICICICIICICI ICICI 


n=1:length(x); 


y(n*M)=x; 
y=[y(M:end) zeros(1,M-1)]; 
on---------------- End of function "up"---------------------------- 
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function noise_power=power_of_noise(x,Tb,Ts,EbNo_ dB) 


ORES O SISO ISSCC GCI ICICI CIR IC 


% Noise Power 

% A function to calculate the noise power associated with a 
% given Eb/No (dB), bit length Tb, and sampling period Ts 
% 

% developed by Pavlos Angelopoulos Nov 2003 

% last modified 11/03 

O62 He 2s ae 2 6 he 2 fe 2 he 2 ee 2 6 2 6 22 fe 2 he 26 6 2 fe 2 fe 22 fe 2 he 2 ae ie 2s fe 2 2 2 2s ie 
signal _power=sum(x.”2)/length(x); 
EbNo=10’(EbNo_dB/10); 

Eb=signal_power*Tb; % energy per bit 
No=Eb/EbNo; 

sigma_No=sqrt(No/(2*Ts)); _% std deviation of noise 
noise_power=sigma_No”2; 


function P=house(x) 


ORES ORS ISI SCC ICICI CCCI ICI ICI IACI CAKE 


% A function to compute the householder matrix P such that 
% transformed_x=P*x is a vector whose elements except the 
% first one are equal to 0. 

% x must be defined as a column vector 

% 

% developed by P.Angelopoulos Nov 2003 

% last modified 11/24 

O63 He 2 6 2 6 2 2 he 2 he 2 he he 2 fe 2 fe 2 ee 2 Ae 2 26 2 2 fe 2 he 2 he Ae 2 2 2 fe 2 ae he ae 2s 2s fe 2 2s 2k 
k=length(x); 

el=[1,zeros(1,k-1)]'; 

v=xtsign(x(1))*norm(x,2)*e1; 

P=eye(k)-2*v*v'/(v'*v); 

o------------------- End of function "house"--------------------- 


function data=spredesp(x,code) 


ORES ORS IS SCISSOR CIGISIC ICICI III ICICI ICICI KK 


" 


% A function to spread (chip) a binary signal x with gold code "code 
% or de-spread the received signal by a replica of the "code" 

% 

% developed by P.Angelopoulos Oct 2003 

% last modified 10/23 

Oo 2 He 2 6 He 6 2 2 he 2 he 2 ee 2 6 2h fe 22 6 2 He 2 6 Ae 2 fe 2 fe 2A 2 6 2 he 2 fe 2 2 fe 2 fe 2 2c he 2 ie 2s fe 2 2s aie 2c 
k=length(x); 

m=length(code); 

repeated_code=repmat(code, | ,ceil(k/m)); 

data=x.*repeated_code; 
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